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SECTION  1 


INTRODUCTION 

Much  has  been  heard  about  the  advantages  of  advanced  composite 
materials  which  accrue  because  of  their  high  specific  stiffness  and 
strength,  excellent  damage  tolerance,  and  superior  fatigue  response 
characteristics,  as  well  as  their  corrosion  resistance  in  hostile 
environments.  Laminated  composites  are  in  increasing  demand,  because 
they  can  actually  be  tailored  by  Che  design  engineer  to  suit  almost 
any  particular  application.  These  laminates  are  fabricated  by  stack¬ 
ing  up  plies  or  laminae  of  unidirectional  fibrous  composites,  with 
each  lamina  oriented  in  different  directions  to  achieve  the  required 
stiffness  and  strength. 

With  the  increased  use  of  laminated  composites,  analytical  tools 
are  needed  for  the  prediction  of  the  laminate  behavior.  The  so-called 
classical  laminated  plate  theory  has  long  been  used  by  many  designers 
to  predict  laminate  response.  Occasionally,  others  have  used  finite 
difference  or  two-dimensional  finite  element  methods  to  predict  the 
laminate  response.  These  methods,  however,  are  limited  in  their 
capabilities.  The  classical  laminated  plate  theory,  for  example, 
neglects  certain  stress  components  which  experimental  observations 
have  shown  to  play  an  important  role  in  failure  of  laminated 
composites.  Experimental  observations  have  also  shown  that  many 
laminated  composites  exhibit  significant  amounts  of  inelastic 
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behavior,  and  that  yielding  of  a  ply  or  lamina  within  a  laminate  does 
not  mean  failure  of  the  laminate.  Approximate  methods  typically 
neglect  this  inelastic  behavior.  Furthermore,  the  effect  of  changes 
in  temperature  and  moisture  on  the  material  properties  is  usually 
neglected  in  these  analyses. 

In  the  present  work,  the  full  three-dimensional  problem  of 
laminated  composites  has  been  addressed.  A  plasticity  model  that 
takes  into  account  the  effect  of  varying  material  properties  due  to 
temperature  and  moisture  changes  has  been  developed  to  describe  the 
inelastic  behavior  of  the  laminate  in  three-dimensional  space.  This 
model  has  been  developed  in  such  a  way  as  to  make  it  possible  to 
incorporate  it  in  a  three-dimensional  finite  element  analysis.  A 
computer  program  has  been  developed  for  the  implementation  of  the 
analysis.  Both  mechanical  and  hygrothermal  loading  conditions  can  be 
handled  by  this  three-dimensional  elastoplastic  finite  element 
program. 

In  Section  2,  an  historical  background  of  the  analysis  of  compo¬ 
site  laminates  is  presented.  The  classical  theory  of  laminated 
plates  and  other  analytical  methods  are  reviewed.  Linear  and 
nonlinear  approximate  analyses  cited  in  the  literature  are  also 
discussed  and  studies  of  thermal  stresses  and  interlaminar  stresses 
are  presented. 

In  Section  3,  the  concepts  of  yield  criterion,  hardening  rule, 
and  flow  rule  are  examined  in  view  of  previous  studies  in  the 
literature.  A  plasticity  model  that  accounts  for  varying  temperature 


and  moisture  content,  and  leads  to  a  stiffness  concept  suitable  for 


use  in  a  three-dimensional  finite  element  analysis,  is  then  developed. 
Finally,  the  finite  element  analysis  and  the  main  features  of  the 
associated  computer  code,  developed  as  the  tool  for  this  work,  are 
discussed.  Example  problems  are  presented  with  results  and  discussion 
in  Section  4,  and  a  general  discussion,  summary,  and  suggestions  for 
future  work  are  included  in  Section  5. 


SECTION  2 


LITERATURE  REVIEW 


2.1  Linear  Analyses 

Tsai  [1]  was  among  the  first  investigators  to  apply  the  classi¬ 
cal  laminated  plate  theory  to  composite  materials.  The  classical 
laminated  plate  theory  (LPT)  that  he  presented  in  1964  is  based  on  the 
work  by  Reissner  and  Stavsky  [2]  in  1961.  Even  earlier  works,  such 
as  that  of  Smith  [3]  published  in  1953,  should  also  be  considered  to 
have  contributed  to  the  development  of  the  laminated  plate  theory. 
Smith  assumed  plane  stress  conditions  in  his  work  on  orthotropic 
plywood.  A  detailed  description  of  the  laminated  plate  theory  is 
given  in  Reference  [4].  The  work  by  Tsai  [1]  is  a  point  stress 
analysis  in  that  it  does  not  take  into  account  boundary  conditions. 
Ashton  and  Whitney  [5]  applied  LPT  to  the  problem  of  laminated  plates 
and  presented  solutions  that  satisfied  the  boundary  conditions  for 
specially  orthotropic  plates,  i.e.,  laminated  plates  composed  of 
orthotropic  layers  such  that  the  orthotropic  axes  of  symmetry  of 
each  layer  coincide  with  the  geometric  axes  of  the  plate.  This  type 
of  laminate  does  not  exhibit  coupling  between  the  bending  moment 
resultants  and  the  twisting  curvatures. 

Although  thermal  loading  is  handled  by  the  laminated  plate 
theory,  to  date  few  attempts  [6]  have  been  made  to  consider  varying 
material  properties  due  to  temperature,  and  none  for  humidity 
changes.  Also  a  significant  limitation  of  LPT  is  that  it  cannot 
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consider  the  out-of-plane  stresses  which,  as  will  be  seen  later  in 
this  chapter,  have  proven  to  be  of  great  importance  in  the  analysis 
of  composite  materials.  Some  computer  codes  based  on  LPT  are  still 
widely  used,  e.g. ,  AC3  [7]  and  SQ5  [8], 

Tsai  [1]  developed  an  analytical  method  to  predict  the  elastic 
constants  of  a  lamina  based  upon  the  properties  of  the  fibers  and  the 
matrix  material.  Then,  using  the  laminated  plate  theory,  he 
theoretically  predicted  the  elastic  constants  of  laminates  made  of 
glass/epoxy  plies.  With  the  aid  of  an  experimental  technique,  Azzi 
and  Tsai  [9]  were  able  to  show  good  agreement  between  the  analytical 
predictions  and  experimental  results  for  glass/epoxy  composites. 
Hill's  yield  condition  [10]  for  orthotropic  materials  under  plane 
stress  conditions  was  employed  to  predict  laminate  failure  in  an 
analytical  study  [11]  of  the  strength  of  transversely  isotropic  lami¬ 
nates.  Inelastic  behavior  was  not  included  in  this  study,  however. 

in  1965,  Tsai  [12]  considered  both  mechanical  and  thermal 
loading  in  his  method  for  modeling  the  load  transfer  from  one  lamina 
as  it  failed  to  the  remaining  unfailed  laminae.  Lamina  failure  was 
predcited  using  Hill's  plane  stress  yield  condition  [10].  After  the 
failure  of  a  lamina,  the  stiffness  matrix  was  degraded  by  setting 
the  appropriate  elastic  constants  to  a  small  fraction  of  their 
original  values.  However,  all  laminae  were  assumed  to  remain 
linearly  elastic  to  failure. 

The  concept  of  a  stepwise  reduction  in  load  carrying  capacity 
after  a  lamina  reached  initial  yielding  according  to  Hill's  yield 
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condition  was  used  by  Chiu  [13].  He  noted  that  failure  in  one 
direction  of  a  certain  lamina  did  not  imply  total  failure  of  this 
particular  lamina,  or  failure  of  the  laminate.  Out-of-plane  stresses 
were  not  included  in  this  study,  however,  and  no  account  was  made 
for  hygrothermal  effects. 

In  1967,  Hoffman  [14]  proposed  a  failure  criterion  similar  to 
Hill's  [10]  yield  condition.  His  theory  contained  linear  terms  and 
hence  could  account  for  differing  tensile  and  compressive  properties. 
Good  agreement  between  theory  and  experiment  for  the  ultimate 
strength  of  laminates  was  shown  for  uniaxial  tensile  and  compressive 
tests;  no  account  was  made  of  the  initial  thermal  stresses  which  are 
induced  during  the  fabrication  of  composite  materials,  due  to  the 
cooldown  from  high  curing  temperatures. 

In  1971,  Tsai  and  Wu  [15]  proposed  a  general  strength  criterion 
for  anisotropic  materials.  They  employed  an  analytical  method  that 
assumed  the  stress-strain  curve  to  be  elastic  to  failure.  Limited 
uniaxial  experiments  showed  good  agreement  with  theoretical  predic¬ 
tions  for  ultimate  strength  of  graphite/epoxy.  This  strength  theory 
and  others  will  be  discussed  in  more  detail  when  yield  surfaces  are 
considered  in  the  next  chapter. 

During  the  fabrication  of  composite  laminates,  residual  thermal 
stresses  are  induced  due  to  cooldown  from  the  cure  temperature.  The 
important  effect  of  these  thermally-induced  stresses  has  long  been 
recognized  in  studies  of  unidirectional  composites,  by  investigators 
such  as  Adams  and  Doner  [16],  Foye  [17],  and  Miller  and  Adams  [18,19]. 
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However,  in  spite  of  this  activity,  in  1974  Hahn  and  Pagano  [6]  noted 
that  no  method  was  yet  available  that  would  account  for  these 
stresses  in  composite  laminates.  They  developed  a  method  based  on 
total  stress-strain-temperature  relations.  Strains  were  assumed  to 
consist  of  thermal  strains  and  mechanical  strains, with  the  latter 
depending  also  on  temperature.  Material  stiffnesses  were  assumed  to 
depend  only  linearly  on  temperature.  By  employing  the  laminated 
plate  theory,  however,  they  carried  along  all  its  limitations 
explained  earlier.  Their  method  did  not  consider  nonlinear  material 
effects. 

2.2  Nonlinear  Analyses 

Throughout  this  work,  the  terms  nonlinear  elastic,  inelastic, 
and  elastoplast ic  will  be  used  to  describe  material  response.  Non¬ 
linear  or  nonlinear  elastic  behavior  implies  that  unloading  from  any 
state  of  stress  is  assumed  to  follow  the  stress-strain  curve,  i.e., 
no  permanent  deformation  occurs.  Inelastic  is  used  to  describe  the 
material  behavior  beyond  its  elastic  limit.  The  term  elastoplast ic 
is  used  to  describe  materials  which  exhibit  plastic  (time-independent) 
deformations  beyond  the  elastic  limit. 

With  the  development  of  high  modulus  fibers  in  the  early  1960's, 
and  hence  the  so-called  advanced  composites,  the  use  of  composites 
increased  sharply.  To  be  able  to  fully  utilize  these  new  materials, 
a  full  and  complete  understanding  of  their  behavior  was  needed. 

Among  many  other  needs,  this  increased  the  demand  for  nonlinear 


analyses. 
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An  early  attempt  was  that  of  Petit  and  Waddoups  [ 20 ]  in  1969. 

The  laminate  behavior  was  established  by  a  piecewise  linear  approach, 
and  incremental  application  of  average  laminate  stresses.  At  each 
increment  of  load,  the  incremental  strains  were  determined  and  added 
to  previous  strains  to  determine  the  total  strains  in  each  lamina. 
These  were  then  transformed  to  the  principal  directions  of  the  lamina. 
Lamina  principal  strains  were  then  referred  back  to  the  respective 
stress-strain  curves,  which  were  composed  of  linear  segments,  and 
updated  stiffnesses  were  determined  for  the  next  load  increment. 
Failure  was  assumed  to  occur  when  a  principal  strain  reached  a 
corresponding  ultimate  value  in  any  lamina.  Inelastic  material 
behavior  was  not  considered,  however.  No  accounting  of  thermal 
effects  was  made,  and  being  a  two-dimensional  analysis,  the  important 
role  of  out-of-plane  stresses  on  failure  mode  could  not  be  handled. 

In  a  similar  analysis  by  Hashin,  et.  al.  [21],  the  Ramber g-Osgood 
equation  [22]  was  used  to  represent  the  lamina  stress-strain  curves, 
this  being  a  three-parameter  curve-fit  equation. 

Cubic  spline  interpolation  curve-fit  functions  were  employed  by 
Sandhu  [23]  to  model  the  lamina  stress-strain  curves.  Incremental 
loading  was  used,  the  material  properties  being  updated  each  incre¬ 
ment  using  equivalent  strains.  Sandhu  proposed  that  it  was  erroneous 
to  determine  updated  material  properties  after  each  loading  increment 
based  on  the  individual  lamina  strains,  as  in  the  work  of  Petit  and 
Waddoups  [20].  He  assumed  that  the  incremental  stresses  could  be 
related  to  the  incremental  strains.  The  ultimate  load  carrying 


capacity  of  the  laminate  was  determined  by  the  plywise  application 
of  a  failure  criterion  which  assumed  that  the  strain  energies  under 
longitudinal,  transverse  and  shear  loadings  were  independent 
parameters.  Better  agreement  between  theory  and  experiment  tor 
various  boron/epoxy  laminates  titan  that  reported  by  Petit  and 
Waddoups  [20]  was  also  indicated.  Sand hu  did  not  consider  hyi'o- 
thermal  loading  in  his  analysis,  and  no  yield  criterion  was  used, 
i.e..  Inelastic  effects  were  neglected. 

Nonlinear  behavior  of  laminated  fibrous  composites  including 
thermal  effects  and  temperature-dependent  material  properties  was 
considered  in  a  study  by  Renieri  and  Herakovich  [24  |.  I'sing  a  two- 
dimensional  finite  element  analysis  and  limiting  the  loading  condi¬ 
tions  to  uniaxial  mechanical  loading,  they  studied  both  nonlinear 
and  thermal  effects.  They  showed  that  these  effects  are  important 
in  the  prediction  of  failure  modes  in  composite  laminates,  and  noted 
that  laminate  failure  is  a  function  of  the  laminate  configuration, 
material,  and  type  of  loading.  Unloading  was  assumed  to  follow 
the  stress-strain  curve.  No  yield  criterion  was  used  in  '.his  two- 
dimensional  analysis,  i.e.,  inelastic  behavior  was  neglected,  and 
failure  was  assumed  when  any  of  the  strains  in  the  principal 
directions  of  the  material  readied  its  ultimate  value.  Tie  analvsis 
was  also  restricted  to  symmetric  laminates. 

2.3  Interlaminar  Stresses 

A  characteristic  of  laminated  romposites  is  that  under  in- plum 
uniaxial  or  biaxial  loading  conditions,  a  laminate  after  develops 
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triaxial  effects,  i.e.,  a  laminate  under  in-plane  loading  only  will 
often,  depending  on  its  lay-up,  develop  out-of-plane  stresses. 

These  out-of-plane  stresses  are  termed  interlaminar  stresses.  Their 
effect  is  most  significant  near  free  edges  where  they  attain  high 
values. 

Plane  stress  conditions  assumed  in  the  classical  laminated  plate 
theory  do  not  allow  for  the  calculation  of  interlaminar  stresses. 
Rather,  only  the  stresses  in  the  plane  of  the  laminate  are  considered. 
The  laminated  plate  theory  is  thus  incapable  of  providing  predictions 
of  some  of  the  stresses  that  actually  cause  failure  of  composite 
laminates,  as  shown  by  experimental  evidence  [25].  In  fact,  inter¬ 
laminar  stresses  are  one  of  the  mechanisms  that  uniquely  characterize 
failure  in  composite  laminates.  These  stresses  may  cause  delamination 
near  a  free  edge,  whether  it  be  at  the  edge  of  a  plate,  around  a 
hole,  etc.  Their  effect  near  an  edge  is  known  as  the  free  edge 
effect . 

In  an  effort  to  evaluate  interlaminar  stresses  and  the  influence 
of  the  stacking  sequence  on  laminate  strength  under  uniform  axial 
loading.  Pipes  and  Pagano  [26]  used  a  finite  difference  technique. 
Their/ study  was  limited  to  angle-ply  laminates,  i.e.,  laminates  with 
an  even  number  of  layers,  with  each  lamina  or  ply  alternately 
oriented  at  +0  and  -6  to  one  of  the  geometric  axes  in  the  plane  of 
the  laminate.  They  noted  that  significant  interlaminar  shear  stresses 
were  required  to  allow  shear  transfer  between  the  plies  of  the  lami¬ 
nate.  In  addition,  the  interlaminar  normal  stress  was  found  to  be  an 
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edge  effect,  restricted  to  a  region  near  the  free  edge  approximately 
equal  to  the  laminate  thickness.  They  reported  strong  evidence  of  a 
singularity  in  interlaminar  normal  stress  at  the  intersection  of  an 
interface  and  the  free  edge.  They  also  suggested  that  such  high 
stresses  in  the  neighborhood  of  the  free  edge  may  be  expected  to 
cause  delamination  of  the  laminate,  in  particular  under  fatigue 
loadings.  Experimental  observations  of  this  phenomenon  had  been 
reported  earlier  by  Foye  and  Baker  [27].  Similar  results  were  also 
reported  by  Pipes  and  Daniel  [28],  in  their  work  on  graphite/epoxy. 
Interlaminar  stress  results  for  cross-ply  graphite /epoxy  laminates, 
i.e.,  laminates  with  alternating  0°  and  90°  plies,  and  other  laminate 
configurations,  subjected  to  uniform  axial  strain  loading  are 
reported  in  Reference  [29]. 

Finite  difference  and  finite  element  methods  [30-42]  both  have 
been  used  to  determine  the  interlaminar  stresses  due  to  in-plane 
loading  of  laminates.  Relatively  little  work  has  been  done  to 
determine  these  stresses  due  to  hygrothermal  loading  or  nonlinear 
effects,  in  spite  of  their  importance. 

Among  the  first  to  use  a  finite  element  analysis  to  determine 
interlaminar  stresses  were  Puppo  and  Evensen  [30],  and  Tsakson  and 
Levy  [31].  Those  investigators  assumed  plane  stress  conditions, 
and  the  fibrous  laminae  were  considered  to  be  orthotropic  layers 
separated  by  isotropic  laminae  of  finite  thickness  that  developed 
only  interlaminar  shear  stresses.  Due  to  the  assumed  plane  stress 
conditions,  however,  interlaminar  normal  stresses  could  not  be 
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modeled.  This  work  and  that  of  Reference  [30]  assumed  linearly 
elastic  material  behavior.  Levy,  et.  al.  [32],  extended  the  work  of 
Isakson  and  Levy  [31]  to  include  plastic  deformation  of  the  shear 
layer,  but  due  to  the  use  of  two-dimensional  analysis  and  the  assumed 
plane  stress  conditions,  they  were  again  unable  to  handle  the  inter¬ 
laminar  normal  stresses.  These  studies  were  also  limited  to 
mechanical  loadings. 

Herakovich  and  Brooks  [33]  considered  laminates  in  the  elastic 
range  under  uniform  axial  strain  loading.  They  used  a  finite  element 
analysis  [27],  and  finite  element  subsections  at  the  free  edge  were 
considered,  i.e.,  the  finite  element  solution  was  reapplied  to  smaller 
sections  near  the  free  edge  to  obtain  more  detail  in  this  region. 
Relatively  high  interlaminar  stresses  at  the  free  edge  and  at  the 
lamina  interface  were  predicted.  They  also  pointed  out  the  signifi¬ 
cance  of  the  interlaminar  stresses  with  respect  to  the  applied  loads 
and  the  stacking  sequence,  and  noted  that  the  magnitude  of  the 
interlaminar  stresses  could  influence  the  strength  of  the  laminate. 
Again  this  analysis  neglected  nonlinear  effects.  To  account  for 
the  interlaminar  normal  stresses,  a  linearly  elastic  three-dimensional 
finite  element  analysis  was  used  by  Rybicki  [34].  He  was  able  to 
consider  both  interlaminar  shear  and  normal  stresses.  The  effects 
of  stacking  sequence  were  studied  under  mechanical  loading  only. 

Very  recently,  Altus,  et.  al.  [35],  presented  a  three-dimensional 
finite  difference  solution  of  the  problem  of  free  edge  effects.  Only 
symmetric  angle-ply  laminates  were  considered  in  this  linearly 
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elastic  analysis.  The  finite  difference  grid  consisted  of  real 
points  inside  the  layer  and  fictitious  points  outside  it  to  satisfy 
the  boundary  conditions  on  the  faces  of  each  layer.  This  study 
confirmed  qualitatively  that  the  normal  interlaminar  stresses 
represent  one  of  the  most  significant  factors  affecting  failure 
of  laminates,  a  result  which  was  also  noted  in  previous  studies 
[36-40].  These  studies  considered  mechanical  loadings  only. 

Results  for  edge  effects  under  mechanical  as  well  as  thermal  loadings 
were  reported  by  Wang  and  Crossman  [41,42].  The  material  properties 
were  assumed  not  to  vary  with  temperature  in  their  study,  however. 

Rybicki  and  Schmueser  [43]  studied  the  effect  of  stacking 
sequence  and  lay-up  angle  on  free  edge  stresses  around  a  hole  in  a 
laminated  plate  under  tension.  A  three-dimensional  analysis  was 
used  to  determine  the  tangential  strain  distributions  around  circular 
holes  in  composite  laminates  under  uniaxial  loading.  Interlaminar 
normal  stress  distributions  around  the  free  edge  of  the  circular 
hole,  changes  in  stacking  sequence,  and  lay-up  angle  were  also 
studied  for  different  graphite/epoxy  systems.  They  did  not  consider 
nonlinear  or  inelastic  material  behavior,  however. 

In  summary,  although  nonlinear  effects  in  advanced  composite 
laminate  systems,  and  both  normal  and  shear  interlaminar  stresses, 
have  been  recognized  to  be  of  major  importance  in  the  prediction  of 
laminate  response  to  mechanical  and  hygrothermal  loading  conditions, 
to  date  no  rigorous,  fully  three-dimensional  nonlinear  analysis  has 
been  developed.  The  aim  of  the  present  work  is  to  do  so. 


SECTION  3 


DEVELOPMENT  OF  ANALYSIS 

Implicit  in  the  development  of  any  plasticity  model  are 
assumptions  associated  with  the  behavior  of  the  material.  Several 
of  these  assumptions  will  be  employed  in  the  present  work;  a  discus¬ 
sion  of  their  implications  follows. 

3.1  Yield  Criterion 

The  initial  yield  condition  defining  the  elastic  limit  of  the 
material  in  a  multiaxial  stress  state  is  referred  to  as  the  yield 
surface  or  the  yield  criterion.  A  discussion  of  the  concept  of 
yield  surface  is  given  in  Reference  [44]. 

The  basic  limitation  of  all  present  yield  criteria  is  the  lack 
of  sufficient  experimental  data  to  support  the  mathematical  model  by 
which  the  yield  surface  is  described .  Phenomenological  effects  such  as 
hardening,  strain  history  dependence,  and  the  Bauschinger  effect  [10] 
are  difficult  to  determine,  and  vary  among  different  materials.  The 
Bauschinger  effect  accounts  for  having  a  different  yield  stress  in 
compression  than  in  tension  after  reversing  the  stress  vector  from 
ary  tensile  plastic  stress  state,  and  vice  versa,  as  seen  in  Figure  1. 
Since  the  exact  role  the  Bauschinger  effect  plays  in  the  behavior  of 
composite  materials  is  still  unknown,  this  effect  is  often  neglected 
in  dealing  with  problems  in  plasticity.  It  will  be  neglected  in  the 
present  analysis  also. 
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Figure  1.  Bauschinger  Effect 

Among  the  most  frequently  used  yield  criteria  are  those 
attributed  to  Tresca  [45]  and  von  Mises  [46].  The  Tresca  condition 
states  that  inelastic  action  at  any  point  in  a  body  begins  only  when 
the  maximum  shearing  stress  on  some  plane  through  the  point  reaches 
a  value  equal  to  the  maximum  shearing  stress  in  a  uniaxial  tension 
specimen  at  yield.  The  more  popular  von  Mises  criterion  [46]  has  been 
extended  to  describe  anisotropic  materials  [10,47-50],  In  the  work 
by  Yamada  [48],  the  yield  condition  was  assumed  to  be  quadratic  in 
the  stress  components,  and  to  reduce  to  the  von  Mises  criterion  when 


the  degree  of  anisotropy  was  small.  The  Hill  yield  criterion  [10] 
took  the  form 


f  (a,  )  =  N.  .o.a  .-1  =  0 
k  i]  i  J 


(1) 


where  N..  are  constants  related  to  the  six  yield  stresses  in  the 
il 

principal  anisotropic  directions. 

Implicit  in  Hill's  criterion  are: 

1)  The  principal  axes  of  anisotropy  do  not  rotate  during 
plastic  flow. 

2)  Orthotropic  material  behavior  is  assumed. 

3)  The  principal  axes  of  anisotropy  either  coincide  with  the 
principal  stress  axes  or  the  transformation  is  known. 

4)  The  anisotropic  parameters  N„  remain  unchanged  during 
plastic  flow. 

Baltov  and  Sawczuk  [51]  generalized  Hill's  yield  criterion  to 
include  combined  isotropic  and  kinematic  hardening.  (Both  will  be 
presented  in  a  later  section) .  They  assumed  a  form 


f(o  )  =  N  (o.  -  a  )(o.  -  a.)  -1=0 
k  ij  i  1  J  3 


(2) 


where  and  a  are  the  kinematic  hardening  parameters.  The  aniso¬ 
tropic  parameters  N^.  are  prescribed  functions  of  the  plastic  strain, 
and  hence  change  during  the  course  of  plastic  flow. 

Several  investigations  have  been  concerned  with  defining  a 


strength  criterion  specifically  for  composite  materials.  A  criterion 
of  this  type  is  used  as  a  failure  surface  or  condition  in  stress 
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space.  When  the  stress  state,  expressed  as  some  function  of  the 
stress  components,  lies  on  this  surface,  it  satisfies  the  failure 
condition.  That  is,  the  material  is  assumed  to  have  failed.  These 
surfaces  could  also  be  used  to  define  yielding,  if  the  parameters 
describing  such  surfaces  were  expressed  in  terms  of  material  yield 
stresses  instead  of  material  ultimate  strengths. 

One  criterion  of  the  type  described  above  is  the  Tsai-Hill 
criterion  [52].  In  fact,  this  strength  criterion  is  based  on  Hill's 
yield  criterion  [10].  Tsai  related  the  anisotropic  parameters  to  the 
usual  failure  strengths  o^,  o^,  and  of  a  lamina.  For  plane  stress 

conditions,  it  takes  the  form 


,  u.  2  .  u.  2 
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Considerable  interaction  exists  between  the  failure  strengths  o^, 
and  of  this  theory  and  those  of  other  criteria  which  assume  that 
axial,  transverse,  and  shear  failures  occur  independently.  A  typical 
element  of  a  fibrous  composite  lamina  Is  shown  in  Figure  2.  The 
three  principal  axes  are:  1  in  the  direction  of  the  fibers, 

2  transversely  in  the  plane,  and  3  normal  to  the  plane  of  the 
lamina. 

Another  strength  criterion  was  presented  by  Tsai  and  Wu  [15]. 
Their  assumption  was 


f(°k)  =  Vi  +  W’j 


1 


(4) 
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Figure  2.  Unidirectional  Fibrous  Lamina 


Although  Ll  seems  ptiysically  attractive,  some  components  in  the 


tensor  polynomial  require  difficult  and  expensive  biaxial  testing  for 
their  determination.  Narayanaswami  and  Adelman  [53]  proposed  a 
modified  criterion,  based  on  the  Tsal-Wu  theory,  in  which  terms 
requiring  biaxial  testing  were  assumed  zero.  This  criterion  could 
also  be  modified  and  used  as  a  yield  criterion. 

3.2  Hardening  Rule 

The  stress-strain  response  after  initial  yielding  differs  among 
various  materials,  and  it  also  differs  among  various  plasticity 
theories.  This  post  yielding  response,  called  hardening,  is 
described  by  specifying  a  subsequent  yield  surface,  often  termed 
the  loading  surface  or  the  hardening  rule.  An  account  of  the  various 
hardening  rules  currently  used  is  given  by  Armen  [54], 

The  choice  of  a  hardening  rule  will  depend  on  its  ability  to 
represent  the  hardening  behavior  of  the  material  under  consideration, 
and  upon  the  ease  with  which  it  can  be  applied.  These  requirements, 
together  with  the  necessity  of  maintaining  mathematical  consistency 
witli  the  yield  function  described  in  the  previous  section,  constitute 
the  criterion  for  the  selection  of  a  hardening  rule. 

Hill  [10]  and  Hodge  [55|  proposed  the  isotropic  hardening  rule, 
which  assumed  that  during  plastic  flow  the  loading  surface  expanded 
uniformly  about  the  origin  in  stress  space,  maintaining  the  same- 
shape,  center  and  orientation  as  the  yield  surface.  Figure  3 
illustrates  the  yield  and  loading  surfaces  when  the  stress  state 
shifts  from  point  1  to  point  2.  Unloading  and  subsequent  reloading 
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in  the  reverse  direction  will  result  in  yielding  at  the  stress  state 
represented  by  point  3.  The  path  2-3  will  be  elastic,  and  0-2  is 
equal  to  0-3.  The  isotropic  representation  of  work  hardening  does 
not  account  for  the  Bauschinger  effect.  Since  fibrous  composites 
are  strongly  anisotropic,  the  directions  of  anisotropy  remain 
effectively  unchanged  during  deformation.  Thus,  the  isotropic  work 
hardening  rule  is  suitable  for  these  materials. 


Figure  3.  Isotropic  Hardening 

The  hardening  behavior  postulated  in  the  theory  of  kinematic 
hardening  is  due  to  Prager  [56,57],  and  assumes  that  during  plastic 
deformation  the  loading  surface  translates  in  stress  space,  main¬ 
taining  the  size,  shape,  and  orientation  of  the  yield  surface.  The 
primary  aim  of  this  theory  is  to  provide  a  means  of  accounting  for 
the  Bauschinger  effect.  Figure  4  is  an  illustration  of  kinematic 
hardening.  The  yield  surface  and  loading  surface  are  shown  in  this 


21 


figure  for  a  shift  of  the  stress  state  from  point  1  to  point  2. 

The  translation  of  the  center  of  the  yield  surface  is  denoted  by  a 


i  j " 


Figure  4.  Kinematic  Hardening 

A  model  of  combined  kinematic  and  isotropic  hardening,  in  which 
the  subsequent  loading  surfaces  expand  and  translate,  has  been 
presented  by  Hodge  [58|.  However,  since  there  are  no  published 
experimental  data  that  show  that  this  model  fits  the  behavior  of 
composites,  its  complexity  does  not  justify  its  use.  Other  hardening 
rules,  proposed  initially  for  metals,  are  discussed  in  References 
[59-631.  To  date,  there  is  no  universally  applicable  model  to 
describe  all  aspects  of  nonlinear  material  behavior.  There  are, 
however,  some  models  better  suited  to  specific  needs  than  others. 

For  a  particular  class  of  problems,  the  preferred  model  is  that  which 
best  combines  mathematical  and  computational  simplicity  with  a  proper 
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representation  of  experimentally  observed  behavior. 


3.3  Flow  Rule 

Plastic  strain  increments  are  related  to  their  corresponding 
stress  increments  by  means  of  a  flow  rule.  A  general  approach  for 
determining  a  flow  rule  is  the  use  of  the  concept  of  a  plastic 
potential.  The  assumption  is  made  that  there  exists  a  scalar 
function  of  stress  f(o^),  from  which  the  components  of  plastic 


strain  increment  are  proportional  to 


Sf 


So 


ij 


If  f(o^)  represents 


the  yield  surface  in  stress  space,  then  the  above  assumption 
represents  a  result  of  Drucker's  postulate  [64],  which  states  that 
the  work  done  by  an  external  agency  during  a  complete  cycle  of 
loading  and  unloading  must  be  non-negative.  Furthermore,  this 
assumption  leads  to  an  incremental  or  associated  linear  flow  theory 
of  plasticity,  in  which  the  increment  of  plastic  strain  is  in  the 
direction  of  the  outward  normal  to  the  surface  represented  by  f(o^) 
in  stress  space,  at  the  current  value  of  stress.  A  strain  rate 
vector  deviating  from  the  outward  normal  to  the  yield  surface  in  a 
direction  independent  of  the  stress  rate  vector  constitutes  the 
nonassociated  linear  flow  theories  of  plasticity.  The  nonassociated 
flow  theories  are  particularly  suitable  for  work-softening  materials 
[65],  and  can  reasonably  fit  the  behavior  observed  in  some  soils. 

An  associated  flow  rule  that  is  generally  used  to  describe 
elastic-plastic  behavior  is  the  Prandtl-Reuss  relation,  which  is  a 
generalization  of  the  Levy-Mises  equations.  The  Prandtl-Reuss 
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assumption  is  that  the  plastic  strain  increments,  dc|j\,  are  propor¬ 
tional  to  the  corresponding  stress  deviator,  o^,  where  the 
instantaneous  non-negative  value  of  the  constant  of  proportionality 
is  dependent  on  the  work  hardening.  The  concept  of  an  effective 
stress  and  effective  plastic  strain  [66]  is  in  itself  an  assumption 
that  is  usually  introduced  to  reduce  the  complexity  of  a  multiaxial 
situation  to  one  that  can  be  related  to  uniaxial  behavior.  Thus, 
the  proportionality  parameter  can  be  the  ratio  of  the  effective 
stress  to  the  effective  plastic  strain  increment. 


3. A  Stiffness  Concept 

The  three-dimensional  problem  of  laminates  of  fibrous  composites 
will  be  formulated  and  analyzed.  The  stiffness  concept  that  leads 
to  the  finite  element  analysis  will  then  be  discussed  in  relation  to 
the  solution  of  elastoplastic.  problems  in  three-dimensional  space. 

A  transversely  isotropic  unidirectional  composite  lamina  will  be 
assumed.  If  the  three  principal  axes  of  anisotropy  are  taken  to  be; 

1  in  the  direction  of  the  fibers,  2  transversely  in  the  plane,  and 
3  normal  to  the  plane  of  the  lamina,  as  shown  in  Figure  2,  the  2-3 
plane  will  be  the  plane  of  transverse  isotropy.  As  discussed  earlier, 
the  directions  of  anisotropy  will  not  change  during  deformation.  A 
quadratic  form  in  the  six  components  of  stress,  similar  to  Hill's 
yield  condition  [10],  can  then  be  chosen  in  the  form 
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2f (o . .)  =  F(a  -  a.)2  +  G(o  -  a  )2  +  H(o.  -  o„)2 
ij  2  3  3  1  1  2 


+  2L  t23  +  2M  t231  +  2N  =  1 


(5) 


where  F,  G,  H,  L,  M,  and  N  are  parameters  characteristic  of  the 
current  state  of  anisotropy. 

In  the  present  study,  the  parameters  of  anisotropy  are  allowed 
to  vary  with  changes  in  temperature  and/or  moisture  content.  It  can 
easily  be  shown  [10J  that 


1 


(a[)2 


G  +  H  ,  2F  = 


1  +  1 


.  y.  2  ,  y.2  ,  y.  2 

(o2)  (op  (Oj) 


1 


(o^)2 


=  H  +  F  ,  2G  = 


(o’)2 


(Oj)2 


(o’)2 


(6) 


1  =  F  +  G  ,  2H  =  — 0  +  1 


(o^)2 


,  y.2  ,  y.2  ,  y.2 

(op  (op  (op 


y  y  y 

where  o',  o',  and  o'  are  the  tensile  yield  stresses  in  the  1,  2,  and 

y  y  y 

3  directions  of  anisotropy.  Also,  if  x'^,  xp  and  x'2  are  the  yield 
stresses  in  shear  with  respect  to  the  principal  axes  of  anisotropy, 
then 


21,  = 


1 


(T^)2 


2M  = 


1 


(t31>2 


2N  = 


(0’2,2 


(7) 


The  form  of  Eq.  (5)  is  valid  only  when  the  principal  axes  of 
anisotropy  are  taken  to  be  the  axes  of  reference;  otherwise  the 
stress  components  must  be  transformed  in  the  manner  described  in 
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Appendix  A.  The  functional  dependence  of  the  parameters  of 
anisotropy  on  temperature  and  moisture  follows  from  Eqs.  (6)  and 
(7)  when  the  yield  stresses  are  expressed  as  functions  of 
temperature  and  moisture  content. 

The  obvious  association,  implied  by  the  term  'work-hardening,' 
between  the  work  used  to  produce  plastic  flow  and  the  hardening 
created,  suggests  the  hypothesis  that  the  degree  of  hardening  is  a 
function  only  of  the  total  plastic  work,  and  is  otherwise  independent 
of  the  strain  path.  The  external  work  dW  per  unit  volume  done  on 


the  element  during  an  infinitesimal  increment  of  strain  de 


the  continued  loading  of  an  element  of  the  material  is 


with 


dW 


a.  .de 
ij 


ij 


A  part  of  this  work 


(8) 


dW 

e 


a,  .de 

ij 


e 

1.1 


(9) 


represents  recoverable  elastic  energy;  the  remainder  is  the  plastic 
work  per  unit  volume,  i.e.. 
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W  =  o. .de*\ 
ij  ij 


(11) 


In  order  for  plastic  work  to  be  performed,  the  state  of  stress 
must  be  on  the  yield  surface,  i.e.,  the  stress  state  must  also 
satisfy  the  condition  given  by  Eq,  (5).  To  enforce  this  constraint, 
the  Lagrange  multiplier  dA  is  used  [67].  Then 


v  —  [o.  .de?.  -  f(o  )dA]  =  0 
3a  i]  i]  ij 


(12) 


which  gives 


de ,  . 


3f 


3o 


dA 


(13) 


1-1 


With  the  use  of  Eq.  (5),  a  set  of  equations  for  the  plastic  strain 
increments  can  then  be  written  as  follows: 


deu 


°ijdX,  i  =  j 


2t  ,  .dX,  i  /  j 
ij  J 


(14) 


where 


=  [H (°1  -  a2)  +  G(o1  -  o3)]/(F  +  G  +  H) 


a2  =  [F(o2  -  o3)  +  H(a2  -  ai)]/(F  +  G  +  H) 


a3  =  [G(o3  -  Oj)  +  F(o3  -  o2)]/(F  +  G  +  H) 


T23  =  Lt23/(F  +  G  +  H> 


(15) 


jr  •»- 
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X31  =  Mt31/(F  +  g  +  h) 


T12  =  Nti2/(F  +  G  +  H) 


Separating  the  strains  into  elastic  and  plastic  components 


gives 


de^  =  +  S^dc^  +  S^da^  +  dc 


de„  =  S21da1  +  s22do2  +  S23d°3  +  d£2 


(16) 


d£3  =  S31d°l  +  S32d°2  +  S33d°3  +  d£: 


and  for  the  shear  components 


dY23  =  S44dT23  +  dY23 


dY31  =  S55dT31  +  dY31 


(17) 


dY12  =  S66dT12  +  dY12 


where  S^.  are  the  coefficients  of  the  elastic  compliance  matrix  [S 


(del  =  [SeWl"  (18) 

The  inverse  of  K<i .  (18)  is 

(do  }  =  [Ce ] {dt  r 
0 

where  [C  ]  is  the  elastic  stiffness  matrix. 


(19) 
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Relating  the  six  parameters  of  anisotropy  to  the  strain  history 

is  an  extremely  complicated  problem.  The  problem  can  be  simplified, 

however,  by  the  assumption  that  the  yield  stresses  must  increase  in 

proportion  with  strain  hardening.  This  assumption  is  justified  by 

the  fact  that  the  directions  of  anisotropy  in  fibrous  composites 

remain  effectively  the  same  during  deformation,  and  that  for  lack  of 

experimental  data,  the  choice  of  a  complex  hardening  rule  is  not 

justified.  If  X  ,  Y  ,  etc.,  are  the  initial  yield  stresses,  then 

according  to  the  assumption  of  isotropic  hardening  above,  X  =  hX^, 

Y  =  hYQ,  etc.,  where  h  is  a  parameter  increasing  monotonically  from 

unity  which  expresses  the  amount  of  hardening.  The  anisotropic 

parameters  must  then  decrease  in  accordance  with  Eq.  (6)  as 
2 

F  =  FQ/h  ,  etc.  The  way  in  which  h  varies  with  the  strain  can  be 
explained  by  analogy  with  the  isotropic  theory  due  to  von  Mises. 

Let 


F  +  G  +  H 


F(o2  -  a3)2  +  G(a3  -  o/  +  H(o1  - 


F  +  G  +  H 


2l4  +  2m4  +  2N,i2 


F  +  G  +  H 


(20) 


be  a  nondimensional  measure  of  the  equivalent  stress  o'.  By  analogy 
with  the  von  Mises  criterion  for  isotropic  materials.  Hill  [10] 
suggested  that  if  there  is  a  functional  relation  between  a  and  the 
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work  W  (this  is  yet  to  be  demonstrated  by  experiment) ,  there  must  be 


one  between  cT  and  the  effective  (or  equivalent)  strain, 
defined  by 


dc ,  as 


de 


■[» 


(F  +  C  +  H) 


1/2rF(Gde9-Hde3)2+G(Hde3-Fde1)2+H(Fde1-Gde2)2 


(FG  +  GH  +  HF)‘ 


2(dy  ) 2  2(dy  )2  2(dy  J*1 
+  - —  +  - — ^ —  + 


M 


N 


1/2 


(21) 


where  de^  and  dy„  are  given  by  Eqs.  (16)  and  (17).  This  is  the 
analogue  of  the  equivalent  stress-equivalent  strain  curve  for 
isotropic  materials,  the  area  under  which  is  equal  to  the  work  per 
unit  volume.  Accordingly, 


dW  =  o(de  -  de  S  )  =  ode  P 

P 


(22) 


But,  from  Eq.  (10) 


dWp  =  +  ...  +  x23dyP3  +  .. 


(23) 


Substituting  for  de^  and  dyP^.  from  Eq.  (14)  into  Eq.  (23)  yields 


2-2 

dW  =  f  a  dX 
P  3 


(24) 


If  an  effective  stress-effective  plastic  strain  curve  is  then 
constructed,  the  slope  of  such  a  curve  at  any  point  will  be 


H’ 


dq 
cfe  P 


(25) 
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from  which 


deP 


do 

H' 


(26) 


Now,  substituting  for  d  "e  p  in  Eq.  (22),  and  equating  the  result  to 

Eq.  (24)  since  both  are  equal  to  the  plastic  work  per  unit  volume 

dW  , 

P 


a  2  _  2 

■tri  do  =  -r  o  dX  =  dW 

ri  J  p 


(27) 


Rearranging 


2  4  2 

j  ado  =  ~  a  H'dX 


(28) 


The  left-hand  side  of  Eq.  (28)  is  recognized  as  the  differential  of 
Eq.  (20),  defining  o'.  Thus, 


3  =  FTTTh  rG(03  '  °l)dal  +  H(01  ~  °2)d0l] 


+[F(a2  -  CTj) ^a2  “  H(a^  -  o2)do2] 


+t-F(a2  -  o3)da3  -  0(a3  -  a^do.^] 


+  2Lt23dT23  +  2MT31dt31  +  2NT12dTi2 


With  the  use  of  the  definitions  of  Eq.  (15),  Eq.  (29)  can  be 


(29) 


rewritten  as 
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2  is  is  is 

-  adF  =  o^do^  +  °2do2  +  °3do3 


+  2T23dT23  +  2T31dT31  +  2T12dT12 


(30) 


To  arrive  at  a  relation  between  stress  and  strain,  Eqs.  (16)  and 
(17)  must  be  inverted.  Rewriting  these  equations  in  matrix  form. 


{de }  =  [Se]{do}  +  {deP> 


Thus, 

[Se]_1{dc>  =  [Se]-1[SG] {do }  +  [Se]-1{deP} 

But  [Se]  ^  =  [Ce].  Also,  because  of  the  symmetry  of  [Se], 

[Se]  1 [ Se ]  =  [I],  where  [I]  is  the  identity  matrix.  Hence,  the  above 
expression  becomes 

[C6] {de  }  =  {do}  +  [C6 ] {dcP } 


or 


{do}  =  [Ce] {de }  -  [Ce  ]  (dcP  } 


Substituting  from  Eq.  (14)  for  d c  , 


*3- 


ij 


'll 


ij 


or 


{dotj}  =  [C^ j 1 (de i  -  d X { A } 


(31) 


l 


(32) 


where,  for  an  orthotropic  material,  i.e.,  a  material  with  three 
planes  of  symmetry. 
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Finally,  substituting  for  dX  from  Eq.  (33)  into  Eq.  (32)  yields  the 
desired  form  for  the  stress-strain  relation 


{da}  =  [CP] (de } 


(37) 


where 


is  the  plastic  stiffness  matrix 
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3.3  Material  Model 

To  apply  the  method  of  analysis  developed  in  the  previous 
section  to  fibrous  composites,  the  material  properties  in  the  1,  2, 

i 

and  3  directions  are  needed.  If  the  material  is  transversely 
isotropic,  the  properties  in  the  2  and  3  directions  are  the  same. 

For  mathematical  consistency  with  the  formulation,  a  relation 
between  the  effective  stress  and  effective  strain  is  required.  Also, 
the  tangent  modulus  H'  of  the  effective  stress-effective  strain 
curve  is  needed,  as  shown  in  Eq.  (36).  Furthermore,  the  dependence 
of  the  material  properties  on  temperature  and  moisture  content  is 
required  if  hygrothermal  loadings  are  to  be  handled,  and  the  actual 
material  response  under  varying  conditions  of  environment  is  to  be 
considered.  Both  the  elastic  and  the  plastic  response  of  a  material 
change  with  temperature,  and  for  polymeric  materials,  with  moisture 
also,  as  evidenced  by  numerous  experiments  [6,16,17].  Therefore, 
the  entire  effective  stress-effective  strain  curve  of  the  material 
for  each  state  of  the  environment  must  be  available. 

Ramberg  and  Osgood  [22],  used  a  three-parameter  model  to 
describe  the  stress-strain  curve  of  a  material.  In  their  model,  the 
strain  was  expressed  as  a  function  of  the  stress,  with  the  three 
parameters  selected  to  best  fit  the  empirical  data.  Richard  and 
Blacklock  [68,69]  developed  another  three-parameter  model  which 
was  found  to  fit  stress-strain  curves  for  composite  materials  more 
accurately  [19].  This  model  is  of  the  form 
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a  =  Ee/[1  +  |Ee/ao|n]1/n  (39) 

where  and  n  are  two  independent  parameters,  and  E  is  the  initial 

slope  of  the  stress-strain  curve.  Since  the  shape  of  an  effective 
stress-effective  strain  curve  is  similar  to  a  tensile  stress-strain 
curve,  a  similar  equation  for  the  effective  stress-effective  strain 
can  be  written  as 


o  =  E  e/[l  +  |e  e/oo|n]1/n  (40) 

where  <j  is  the  effective  stress  and  ~e  is  the  effective  strain,  as 
defined  by  Eqs.  (20)  and  (21),  respectively.  The  two  independent 
parameters  <Jq  and  n,  together  with  the  third  parameter  E,  which  is 
the  initial  slope  of  the  curve,  are  selected  to  best  fit  the 
experimental  data. 

The  tangent  modulus  H'  of  the  effective  stress-effective  plastic 
strain  curve  is  related  to  E  as  [70] 


(41) 


where  E  is  found  by  differentiating  Eq.  (40)  with  respect  to  r. .  The 
resulting  equation  is 


1+n 

et  -  E/[i  +  |e  r/oQ|ni  n 


(42) 
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By  fitting  Eq.  (40)  to  the  effective  stress-effective  strain 
curves  for  different  temperatures  and  moisture  contents,  a  functional 
relationship  of  the  parameters  E,  F  ,  and  n  to  temperature  and 
moisture  can  be  established.  In  a  similar  manner,  functional  rela¬ 
tionships  can  also  be  found  for  all  other  material  properties. 

3.6  Finite  Element  Formulation 

The  finite  element  technique  is  now  widely  accepted  as  an 
accurate  method  of  stress  analysis.  Progress  has  been  made  on 
three  fronts,  all  of  which  contribute  to  the  strength  and 
flexibility  of  the  method.  First  of  all,  the  relation  of  the  finite 
element  method  to  previous  well-established  methods  in  continuum 
mechanics  has  given  it  a  firm  foundation.  Secondly,  the  search  for 
and  development  of  the  many  types  of  finite  elements  that  suit 
continuity  and/or  equilibrium  requirements  has  given  it  a  wide  area 
of  application.  Finally,  extension  of  the  method  to  the  study  of 
both  material  and  geometric  nonlinearities  has  resulted  in  more 
realistic  models  and  design  methods. 

Prior  to  the  widespread  application  of  digital  computers,  the 
inelastic  behavior  of  solids  was  one  of  the  most  intractable 
problems  in  solid  mechanics.  The  combination  of  large  digital 
computers  and  the  finite  element  method  has  made  possible  the 
solution  of  most  of  these  hitherto  intractable  problems.  However, 
in  most  of  the  previous  literature  on  the  subject,  simple  constant 
stress  elements  have  been  used.  This  was  largely  motivated  by  the 
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difficulty  of  performing  an  exact  integration  in  an  element  contain¬ 
ing  both  elastic  and  plastic  regimes.  With  numerical  integration 
techniques  already  introduced  into  the  formulation  of  complex 
elements  [71-73]  in  the  analysis,  such  difficulties  disappear.  Among 
the  different  numerical  integration  techniques,  the  Cause  method 
[70]  is  most  widely  used. 

A  recent  comparison  of  various  complex  elements  in  three- 
dimensional  analyses  by  Clough  [74]  has  shown  that  isoparametric, 
numerically  integrated  brick  elements  are  more  efficient  in  elasto- 
plastic  analyses  than  simple  elements.  An  isoparametric  formulation 
in  three-dimensional  space  is  presented  in  detail  in  Appendix  B. 

This  type  of  element,  viz,  the  isoparametric  element,  will  be  used 
in  the  present  work. 

The  data  required  to  incorporate  any  element  into  a  static 
analysis  can  be  conveniently  organized  into  four  characteristic 
matrices.  These  matrices  define  the  elastic  or  plastic  behavior, 
the  spatial  assembly  into  the  structure,  and  the  required  output 
information  for  each  element.  The  four  matrices  are: 

1)  The  element  stiffness  matrix,  [Ce]  or  [CP]. 

2)  The  element  initial  stress-free  IF!  or  {dl  vector. 

mo  mo 

3)  The  element  assembly  matrix. 

4)  The  element  output  matrix. 

The  element  initial  stress-free  vector  depends  on  whether  a 
stress  field  or  a  displacement  field  is  assumed  for  the  element. 


For  any  type  of  element,  a  force-deformation  relationship  is 
needed,  which  leads  to  the  first  two  characteristic  matrices,  i.e. 


the  element  stiffness  matrix  and  the  initial  stress-free  vector. 
The  form  of  this  relation  is 


{F}  =  [C]  {d}  -  {F> 

m  mm  mo 


(43) 


where 

{F}  =  independent  generalized  forces  for  element  m;  a  second 

subscript  zero  denotes  initial  values. 

[C]  =  element  stiffness  matrix  (non-singular), 

m 

{d  }  =  independent  generalized  deformations  for  element  m. 
m 


The  third  characteristic  matrix  is  the  element  assembly  matrix. 
To  assemble  the  structural  stiffness  matrix  in  the  assumed  displace¬ 
ment  method,  which  is  used  in  this  work,  the  independent  deformation 

variables  {d}  for  each  element  have  to  be  transformed  into 
m 

equivalent  nodal  displacements  (A)  in  the  global  system,  using  the 

m 

transformation  (assembly)  matrix  [a]  ,  i.e., 

m 


(d)  =  [a]  (A) 

in  mm 


(44) 


Depending  on  whether  the  assumed  displacement  method  or  the 
assumed  stress  method  is  used,  the  for.ith  characteristic  matrix 
relates  the  independent  variables  to  the  required  output  information, 
i.e.,  stresses  in  the  former  case,  and  displacements  in  the  latter. 
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3.7  Computer  Program 

While  the  principles  of  three-dimensional  elastic  stress 
analysis  by  the  finite  element  method  have  been  obvious  from  the 
early  days  of  the  development  of  the  method  [75-77],  their  practical 
implementation  leads  to  some  immediate  obstacles.  As  the  number  of 
elements  increases,  so  do  the  number  of  degrees  of  freedom,  increasing 
the  size  of  the  stiffness  matrix,  which  requires  larger  computer 
in-core  storage.  In  early  approaches  to  three-dimensional  analyses, 
the  simple  tetrahedral  element  was  the  obvious  choice  [76,78].  The 
adequacy  of  results  of  the  two-dimensional  simple  triangle  gave 
similar  confidence  in  the  tetrahedral  element  [79].  However,  it 
was  soon  realized  that,  although  convergence  to  the  exact  solution  is 
guaranteed  as  the  number  of  elements  is  increased,  it  was  too  slow 
for  problems  of  even  moderate  size  [80].  While  the  tetrahedron  has 
certain  advantages  in  its  formulation  [77],  it  is  nevertheless  an 
inconvenient  shape  to  deal  with  in  grid  generating,  i.e., 
topologically,  and  usually  several  have  to  be  combined  to  form  easily 
managed,  hexahedral  shapes  [79].  Various  families  of  isoparametric 
elements  were  introduced  by  Zienkiewicz,  et .  al.,  in  1967  [7V]. 

These  elements  are  more  efficient  than  tetrahedrons,  as  mentioned 
earlier,  and  will  be  utilized  in  the  present  analysis,  thus  allowing 
a  greater  accuracy  to  be  achieved  for  a  given  number  of  degrees  of 
freedom  and  given  computation  time. 

In  the  present  work,  a  computer  program  has  been  developed 
making  use  of  many  of  the  more  recent  developments  in  both  areas  of 


finite  element  analysis  and  computer  programming.  The  following  is 
a  description  of  the  main  features  of  this  program  'NCLAP' 

(Nonlinear  Composite  Laminate  Analysis  Program) . 

A  modular  approach  has  been  adopted  for  NCLAP,  with  the  various 
main  finite  element  operations  being  performed  by  separate  sub¬ 
routines.  Figure  5  shows  the  organization  of  the  program.  The 
basic  finite  element  steps  are  performed  by  primary  subroutines, 
which  rely  on  auxiliary  subroutines  to  carry  out  secondary 
operations.  The  construction  falls  into  three  phases. 

Phase  1.  Data  are  input  and  checked  for  possible  preparation 
errors,  an  important  feature  when  considering  the  amount  of  data 
input  required  for  three-dimensional  problems. 

Phase  2 .  Stiffness  and  stress  matrices  and  the  applied  load 
vector  are  generated.  The  nature  of  laminated  composite  problems 
requires  elements  with  large  aspect  ratios,  i.e.,  the  ratio  between 
minimum  and  maximum  characteristic  dimensions  of  an  element  in  the 
mesh.  This  is  because  the  thickness  of  a  lamina  is  typically  very 
small  compared  to  its  surface  area.  Large  aspect  ratios  lead  to  an 
ill-conditioned  stiffness  matrix,  i.e.,  the  diagonal  terms  become 
very  small  compared  to  the  off-diagonal  terms,  a  situation  that  can 
lead  to  large  solution  errors.  However,  by  using  reduced  integration 
techniques  [81 J,  the  problem  can  be  overcome. 

Phase  3.  A  'frontal'  solution  technique  [82,83]  is  used  for 
the  solution  of  the  stiffness  equations.  The  advantages  of  using 
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this  method  rather  than  the  well-known  banded  matrix  method  are: 

1)  In  solving  the  stiffness  equations  using  a  banded  matrix, 
the  order  in  which  the  nodes  are  numbered  is  very  important 
since  it  influences  the  bandwidth.  Using  the  frontal 
solution  technique,  however,  does  not  require  any  ordering 
of  nodal  numbers.  Hence,  if  a  mesh  is  to  be  modified  at  a 
later  time,  no  renumbering  is  needed.  This  saves 
considerable  time  and  effort  in  data  preparation. 

2)  For  higher  order  elements,  less  core  storage  is  needed. 
Several  examples  which  justify  this  statement  can  be  found 
in  Reference  [84 ] . 

3)  Since  variables  are  eliminated  in  this  method  as  soon  as 
conceivably  possible,  operations  on  zero  coefficients  are 
minimized  and  the  total  number  of  arithmetic  operations  is 
less  than  with  other  methods.  Thus,  less  storage  and 
computer  time  is  used. 

4)  Because  any  new  equation  occupies  the  first  available  space 
in  the  front,  there  is  no  need  for  a  bodily  shifting  of  the 
in-core  equations  as  in  many  other  large  capacity  equation 
solvers. 

In  any  incremental  analysis,  the  use  of  smaller  load  increments 
implies  a  larger  number  of  increments  to  achieve  the  same  total 
applied  load.  Hence,  more  time  is  spent  in  reconstructing  stiffness 
equations.  A  main  feature  of  the  computer  program  NCLAP  developed 
in  the  present  study  is  that  it  allows  for  the  use  of  small  load 


A3 

increments  without  increasing  the  computing  time  significantly.  The 
stiffness  matrix  is  reconstructed  only  for  those  elements  that 
become  plastic,  or  when  hygrothermal  loading  is  considered. 

A  description  of  the  input  data  required  is  given  in  Appendix  D. 


SECTION  4 


EXAMPLE  PROBLEMS 

The  method  of  analysis  presented  in  this  work,  together  with 
the  computer  program  NCLAP  developed  for  its  implementation,  can  be 
applied  to  a  very  wide  range  of  problems.  This  section  illustrates 
a  few  of  these  possible  applications.  Four  different  problems  will 
be  solved  using  the  analysis  and,  whenever  possible,  comparisons 
with  results  obtained  using  other  methods  will  be  made.  These 
four  problems  cover  the  areas  of  generally  orthotropic  laminated 
beams  and  plates  in  bending,  free  edge  effects  in  laminated  plates, 
and  the  problems  associated  with  these  edge  effects  around  circular 
holes.  Both  mechanical  and  hygrothermal  loadings  will  be  considered, 
and  nonlinear  material  effects  will  be  included. 

4.1  Flexure  of  a  Composite  Beam  Under  Three-Point  Loading 

The  problem  of  a  composite  beam  under  three-point  loading  is 
presented  here  and  results  are  compared  with  results  for  the  same 
problem  using  classical  laminated  plate  theory  as  reported  by 
Adams  and  Miller  [85].  The  composite  beam  is  similar  to  an 
unnotched  standard  Charpy  impact  specimen,  being  simply  supported, 
with  a  span  s  =  32.5  mm,  and  width  b  equal  the  height  h  =  6.5  mm, 
as  shown  in  Figure  6.  This  beam  consists  of  thirteen  plies,  with 
two  plies  at  0°,  plies  at  +45°  and  -45°,  two  0°  plies,  one  ply  at 
90®,  two  plies  at  0®,  plies  at  -45®  and  +45°,  and  two  0°  plies. 
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In  accordance  with  the  standard  notation  of  Reference  [7],  the  above 
composite  laminate  is  identified  as 

[0  /+45/0„/90] 
z  —  z  s 

where  the  bar  over  the  90°  ply  indicates  that  the  plane  of  symmetry 
passes  through  the  middle  of  this  ply.  Symmetric  laminates,  such 
as  the  one  presented,  contain  pairs  of  plies  aligned  at  the  same 
angle  at  equal  distances  on  opposite  sides  of  the  midplane.  For 
this  type  of  laminate  there  is  coupling  between  twisting!  curvatures 
and  bending  moments,  i.e.,  the  D1&  and  D2&  terms  in  the  [D]  matrix 
of  the  classical  laminated  plate  theory  [4]  are  not  zero. 


Figure  6.  Simply  Supported  Beam  Under  Three-point  Loading. 
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The  material  used  was  a  Modulite  5206  graphite/epoxy  composite 
[86],  the  mechanical  properties  of  which  are  as  follows: 


E 


E 
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11 

22 

12 

12 


=  141GN/m  (20.52  Msi) 


E33  =  15  GN/m  (2.15  Msi) 


G23  =  G31  =5-5GN/m  (0*8  Msl) 


V23  V31  °‘25 


(45) 


The  beam  lay-up  is  shown  in  Figure  7,  with  the  z-scale 
expanded  to  show  detail.  The  three-dimensional  finite  element  grid 
used  consisted  of  ten  layers  of  elements  with  six  elements  in  the 
x-direction  and  one  element  in  the  y-direction  in  each  layer.  In 
this  grid,  the  two  0°  plies  were  treated  as  one  layer  and  the  90° 
ply  as  two  layers.  The  90°  ply  was  represented  by  two  layers  of 
elements  so  that  the  deflection  at  the  midplane  could  be  calculated 
and  compared  to  that  of  Reference  [85],  since  only  midplane  deflec¬ 
tions  could  be  calculated  by  the  latter. 

The  load  was  applied  as  concentrated  forces  at  the  nodal  points 
of  the  midspan,  i.e.,  at  x  *  s/2  as  shown  in  Figure  7.  A  total 
force  of  100  N  (22.5  lb)  was  applied  in  the  negative  z-direction  and 
the  deflections  at  x  *  s/2,  at  each  ply  interface,  were  calculated. 
Table  1  gives  the  deflections  at  midspan  for  different  values  of  z. 
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In  the  study  by  Adams  and  Miller  [85],  the  classical  laminated 
plate  theory  was  used  to  solve  the  same  problem,  viz,  bending  of  a 
composite  beam  under  three-point  loading  (Figure  6) .  They  assumed 
a  double  sine  series  for  the  load  in  the  form 
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-  sin  - — — 

nm  s 
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where  p^  is  the  load  intensity,  a  is  the  location  of  the  center  of 
application  of  the  load,  and  2b  is  the  length  of  the  beam  over  which 
the  load  is  applied.  The  value  of  0.02s  was  used  for  b  to  approximate 
a  concentrated  load  at  midspan,  as  illustrated  in  Figure  6.  They 
obtained  a  solution  for  the  deflection  w  in  the  form 


P  s 
V  m 

w  =  l 
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2  2'  2 
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s  .  .  mux 

TT-)'1"  T 

m  tt 


in  which  A,.,,  and  are  terms  of  the  well-known  [A]  and  [ I) ] 
of  the  classical  laminated  plate  theory  [4],  and  is  an 
experimentally  determined  shear  correction  factor.  For  the 
of  this  example,  these  values  are  [85] 


(47) 


mat  r 1 ces 


laminate 


A55  =  36.8  MNm-1 (2.1  x  105  lb  in-1) 

DX1  =  3.2  kN-m  (2.8  x  10A  in-lb) 

k.  s 
Ki 


0.66 
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Using  the  method  of  Reference  {85],  as  described  above,  only 
the  deflection  at  the  midplane  (z  =  h/2)  can  be  calculated;  the 
maximum  deflection  at  midspan,  as  calculated  from  Eq.  (47),  is 

w  =  -7.760pm  (-0.306  x  10_6  in)  (48) 

max 

This  value  corresponds  to  the  value  at  z  =  3.25  mm  in  Table  1,  which 
is 

wi=h/2  =  _8*034hm  (-0.316  x  10_6  in)  (49) 

The  theory  of  Reference  [85],  as  represented  by  the  value  given 

in  Eq.  (48),  was  shown  to  predict  values  of  deflection  somewhat 

lower  than  obtained  experimentally.  Thus,  the  present  analysis, 

represented  by  Eq .  (49),  is  actually  in  better  agreement  with 

available  experimental  data,  possibly  due  to  the  fact  that  the 

and  D„,  terms  were  neglected  in  the  analysis  of  Reference  [85]. 

Table  1  also  shows  the  effect  of  interlaminar  stresses  on  the 

deflection  of  the  midspan.  The  midspan  deflection  is  not  the  same 

at  each  ply  interface;  rather  there  is  a  change  in  the  thickness 

of  each  lamina  due  to  the  concentrated  applied  load  at  midspan. 

This  effect  cannot  be  detected  by  the  classical  laminated  plate 

theory,  and  is  not  predicted  by  Reference  [85).  The  magnitude  of 

the  maximum  interlaminar  normal  stress  o  is  8  percent  of  the  maxi- 

z 

mum  in-plane  stress  o  for  this  example.  This  is  a  significant 
value,  since  the  transverse  stiffness  is  only  4  percent  of  the 


Table  1-  Midspan  Deflections  of  a  Simply  Supported 
f°2^— ‘ ^5/02/90 ]g Laminated  Composite  Beam 

Under  Three-Point  Loading  as  Predicted  by 
the  3-D  Finite  Element  Analysis  (NCLAP) 


Distance  from 
Bottom  Surface 
_z  (mm) 


Deflection 
w  (-pm) 

7.726 


7.797 


7.843 


7.986 


8.034 


51 


longitudinal  stiffness-  The  effects  of  interlaminar  stresses  will  be 
discussed  in  more  detail  in  the  last  two  examples. 

4.2  Cylindrical  Bending  of  Laminated  Plates 

In  the  past,  laminated  plates  have  often  been  analyzed  using  the 

assumption  that  they  behaved  as  specially  orthotropic  laminates. 

For  this  type  of  laminate,  the  D, ,  and  D  terms  of  the  [D]  matrix 

lb  zb 

are  zero.  Most  practical  laminates,  however,  do  not  fulfill  these 
conditions  of  special  orthotropy. 

Unlike  other  methods  in  the  literature,  the  method  of  analysis 
presented  in  the  present  work  car.  handle  any  lay-up  and  is  not 
restricted  to  symmetric  or  specially  orthotropic  laminates.  To  be 
able  to  check  the  results  of  the  present  computer  program,  however, 
problems  reported  in  the  literature  are  modeled  and  results  compared. 
One  such  problem  is  that  of  cylindrical  bending  of  an  antisymmetric, 
specially  orthotropic  (cross-ply)  laminated  plate. 

Shown  in  Figure  8  is  the  three-dimensional  grid  used  to  model 
an  antisymmetric  laminated  plate  of  the  configuration  [0/90/0/90].^. 
The  plate  Is  30  mm  long,  20  mm  in  width,  and  consists  of  four  plies 
each  of  which  is  0.4  mm  thick.  Each  ply  is  represented  by  a  layer 
of  24  elements,  with  6  elements  in  the  x-direction,  and  4  in  the 
y-direction.  Thus,  the  grid  consists  of  a  total  of  9 b  elements, 

175  nodes  and  525  degrees  of  freedom. 

The  plate  is  simply  supported  along  x  =  0  and  x  =  a  and  is  free 
along  the  other  two  edges  (see  Figure  8).  A  constant  transverse 


53 


load  is  applied  by  means  of  nodal  forces.  Each  node  on  the  top 
surface  is  loaded  by  a  IN  (0.225  lb)  force  in  the  negative  z- 
direction,  corresponding  to  a  uniform  transverse  applied  load  of 
q ^  =  0.16  MPa  (23.2  psi) .  The  maximum  deflection  at  the  center  of 
the  plate  as  calculated  by  NCLAP  is 


w|  =  -0.0492  mm  (-1.9370  x  103  in)  (50) 

z=2h 

The  same  problem  was  also  solved  using  the  classical  laminated 
plate  theory  solution  presented  in  Reference  [5],  which  is  based  on 
the  assumptions  that  the  plate  is  of  infinite  length  in  the  y- 
direction  and  that  the  transverse  load  is  independent  of  y,  i.e., 
q  =  q (x) .  Thus,  the  deformed  surface  will  be  cylindrical,  i.e., 


uq  =  uq(x),  vq  =  0,  w  =  w(x) 


(51) 


where  u^  and  v^  are  the  midplane  displacements  in  the  x  and  y 
directions,  respectively.  According  to  Reference  [5],  the  maximum 
deflection  occurs  at  the  center  of  the  plate  and  is  given  by 


w 

max 


IAnV 

384  D 


where 
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(52) 


(53) 


and  w  does  not  vary  with  z,  as  indicated  by  Eq.  (51). 
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With  the  aid  of  the  computer  program  AC3,  which  is  a  point 
stress  analysis  program  [7],  the  [A],  [B],  and  [D]  matrices  of  the 
classical  laminated  plate  theory  were  evaluated  for  the  [0/90/0/90]^ 
configuration.  Using  Eqs.  (52)  and  (53),  the  maximum  deflection  was 
then  calculated  as 

w  =  -0.0653  mm  (2.5708  x  10  ^  in)  (54) 

max 

As  calculated  by  the  present  program,  displacements  in  the  y-direction, 

i.e.,  v  ,  satisfied  Eq.  (51).  Displacements  in  the  x-direction, 
o 

however,  showed  dependence  on  the  other  two  coordinates,  i.e., 
u  =  u(x,y,z).  At  this  point,  the  difference  between  the  value  of 
w  in  Eqs.  (50)  and  (54)  cannot  be  adequately  explained. 

The  w  deflection  variation  along  x  =  a/2  is  given  in  Table  2, 
at  every  ply  interface.  A  slight  decrease  in  w  occurs  away  from 
the  boundaries  y  =  0  and  y  =  b.  This  change  of  w  at  the  boundaries 
is  caused  by  the  edge  effects  due  to  the  presence  of  interlaminar 
stresses.  It  does  not  exist  for  laminates  witli  no  interlaminar 
stresses,  e.g.,  isotropic  laminates.  The  significance  of  inter¬ 
laminar  stresses  will  be  discussed  in  the  following  examples. 

4.3  Interlaminar  Stresses  in  Laminated  Plates 

As  stated  earlier,  laminated  composites  develop  out-of-plane  or 
interlaminar  stresses  even  under  uniaxial  loading.  It  was  also 
shown  that  these  stresses  cause  delamination  of  the  laminate.  Two- 
dimensional  analyses  presented  in  the  literature  cannot  accurately 
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Table  2.  Plate  Deflections  Under  Uniform  Transverse 
Load  Along  x  =  a/2,  as  Predicted  by  the 
3-D  Finite  Element  Analysis  (NCLAP) 


Distance  from 
Bottom  Surface 
z  (mm) 

Deflection  w  (nan)  at  Values 
y/b  Indicated 

of 

0.0 

0.25 

0.50 

0.75 

1.0 

0.0 

0.0492 

0.0489 

0.0488 

0.0489 

0.0492 

0.4 

0.0492 

0.0489 

0.0489 

0.0489 

0.0492 

0.8 

0.0492 

0.0489 

0.0489 

0.0489 

0.0492 

1.2 

0.0492 

0.0489 

0.0489 

0.0489 

0.0492 

1.6 

0.0492 

0.0489 

0.0488 

0.0489 

0.0492 
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handle  the  problem  of  interlaminar  stresses,  since  this  problem  is 
three-dimensional  in  nature.  Results  obtained  with  these  methods 
will  be  compared  to  results  using  the  present  analysis  in  the 
following  example. 

Four-ply  symmetric  laminates  will  be  considered.  Because  of 
symmetry,  only  one  quarter  of  the  upper  two  plies  of  the  plate  need 
be  considered,  as  shown  in  Figure  9.  The  origin  is  therefore  at  the 
midplane.  The  three-dimensional  finite  element  grid  used  consisted 
of  two  layers  of  elements.  Each  layer,  being  30  mm  in  the  x- 
direction  and  20  mm  in  the  y-direction  (see  Figure  9) ,  was  con¬ 
structed  by  six  elements  in  the  x-direction  and  four  in  the 
y-direction.  Ply  thickness  was  taken  to  be  0.4  mm;  thus  a  b/h 
ratio  of  50  was  achieved,  where  b  is  the  width  of  the  plate  in 
the  y-direction  and  h  is  the  lamina  thickness. 

To  show  an  important  and  useful  feature  of  the  present 
analysis,  viz,  the  ability  to  handle  hygrothermal  as  well  as 
mechanical  loadings,  or  even  a  combination  of  these  two  types  of 
loadings,  the  material  properties  over  a  range  of  temperature  and 
moisture  levels  must  be  known.  Few  data  have  been  published 
that  describe  the  dependence  of  material  properties  on  environmental 
changes.  One  such  set  of  data  available  is  for  the  Hercules 
AS/3501-6  unidirectional  graphite/epoxy  composite,  which  will  be 
used  in  this  example.  The  functional  dependence  of  the  properties 


of  this  material  on  temperature  and  moisture  were  experimentally 


determined  at  the  University  of  Wyoming  [87].  These  properties  are 
represented  by  the  second  order  polynomial 

P(T,M)  =  C,T2  +  C„M2  +  C.T'M  +  C.T  +  CcM  +  (55) 

1  J.  J  4  3  0 

where  P  is  the  functional  material  property  of  interest,  T  is 
temperature  (°C) ,  M  is  moisture  content  (weight  percent) ,  and  the 
C's  are  coefficients  given  in  Table  3. 

Although  the  present  method  can  handle  both  applied  nodal 
forces  and  displacements,  most  analyses  presented  in  the  literature 
can  only  handle  uniaxial  strain.  For  this  reason,  i.e.,  to  be 
able  to  compare  results,  the  same  type  of  loading,  viz,  uniaxial 
strain  was  used  here.  For  each  laminate  considered,  a  uniform 
displacement  was  applied  in  the  x-direction  such  that  T  =  0.01 
percent.  The  average  in-plane  stress  cf^  is  also  given  for  each 
case.  This  will  permit  the  comparison  of  the  relative  magnitudes 
of  the  applied  in-plane  stresses  to  the  interlaminar  stresses , to 
ascertain  the  importance  of  the  latter. 

4.3.1  Cross-Ply  Laminates 

The  interlaminar  stress  distributions  for  [90/0]  and  [0/90] 

s  s 

laminates  are  shown  in  Figure  10.  The  interlaminar  normal  stress 

a  ,  plotted  in  Figure  10a,  is  very  small  at  y/b  =  0.  As  y/b 
z 

increases  toward  the  free  edge  of  the  plate,  o increases  (in 
absolute  sense)  to  a  peak  value,  then  reverses  direction.  This 
peak  value  occurs  at  y/b  =  0.5  for  the  [90/0]g  laminate,  and  at 


Table  3.  Coefficients  of  Polynomials  in  Eq.  (55)  for  the  Mechanical 
Properties  of  Hercules  AS/3501-6  Unidirectional  Graphite/ 
Epoxy  Composite  [87] 


Material 

Property 

C 

C 

c. 

c 

C 

c. 

(MPa) 

12  i  4  5  6 

E 

-0.424 

-1.357 

-1.329 

-13.66 

8.781 

1.301  (10J) 

*11 

,  ~4 

E 

-0.164 

-0.211 

-0.133 

-10.97 

0. 121 

1.119  (10  ) 

*22 

4 

E  - 

-0.164 

-0.211 

-0.133 

-10.97 

0.121 

1.119  (10  ) 

*33 

-4 

,3 

-0.046 

-0.095 

1.33  (10  ) 

-5.451 

-4.035 

3.74S  (10  ) 

^23 

-4 

,  3 

-0.086 

-0.178 

2.50  (10  ) 

-10.23 

-7.570 

7.031  (10  ) 

31 

-4 

G  „ 

-0.086 

-0.178 
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Table  A.  Poisson's  Ratio,  and  Thermal  and  Moisture  Expansion 
Coefficients  for  Hercules  AS/3501-6  Unidirectional 
Graphite/Epoxy  Composite  [87] 
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y/b  =  0.75  for  the  I 0/90 ] ^  configuration,  i.e.,  closer  to  the  free 

edge  in  the  latter  configuration.  The  peak  value  of  o for  the 

[0/90]  laminate  is  almost  double  that  reached  in  the  [90/0] 
s  s 

laminate.  At  y/b  =  0.96,  the  interlaminar  normal  stress  shows  a 

high  but  finite  value;  tension  in  the  [0/90]^  configuration  and 

compression  in  the  [90/0]^  laminate.  Thus,  the  distribution  of  the 

interlaminar  normal  stress  in  the  [0/90]  laminate  is  not  a  mirror 

image  of  the  distribution  in  the  [90/0]  laminate.  This  behavior 

s 

agrees  with  results  recently  reported  by  Spilker  and  Chou  [88]  in 
their  work  on  interlaminar  stresses.  A  slightly  different 
graphite /epoxy  material  was  used  in  their  work,  e.g.,  the  longitudinal 
modulus  =  1.4  x  10^  MPa  as  compared  to  a  value  of  1.2  x  105  MPa 
for  the  Hercules  AS/ 3501-6  material  used  in  this  example.  The 
purpose  here,  however,  is  to  show  the  trend  in  the  material  behavior 
rather  than  a  point-by-point  stress  comparison.  Spilker  and  Chou 
also  reported  tha ;  satisfying  the  traction-free  edge  condition,  i.e., 
the  boundary  conditions  at  the  free  boundaries,  would  result  in  the 
convergence  of  all  stresses  to  finite  values  at  the  free  edge,  in 
direct  conflict  with  the  idea  of  stress  singularities  [39,40],  To 
satisfy  the  traction-free  edge  condition,  special  elements  must  be 
used,  such  as  those  of  Reference  [88],  or  the  hybrid  elements  of 
Reference  [79].  These  are  complex  elements,  however,  and  Reference 
[88]  shows  that  satisfying  the  traction-free  edge  conditions,  while 
giving  more  accurate  point-by-point  stress  values,  does  not  change 


the  stress  gradient. 
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As  shown  in  Figure  10b,  the  interlaminar  shear  stress  distri¬ 
bution  t  for  the  [0/90]  laminate  is  also  not  a  mirror  image  of 
yz  s 

that  for  the  [90/0]  laminate.  The  interlaminar  shear  stress  t 

s  yz 

increases  in  absolute  value  from  y/b  =  0  to  y/b  =  0.96,  where  it 
reaches  a  peak  value  near  the  free  edge,  then  drops  to  zero  at  the 
free  boundary  (y/b  =  1.0),  as  dictated  by  the  traction-free  boundary 
condition. 

Figure  11  is  a  plot  of  the  interlaminar  normal  stress  a 

z 

through  the  thickness.  To  increase  information  about  the  distribu¬ 
tion  through  the  thickness,  the  finite  element  grid  was  changed  for 
the  plot  of  Figure  11.  Each  ply  was  represented  by  two  layers  of 
elements  instead  of  one.  It  is  again  noted  that  the  stress 

distribution  for  the  [90/0]  laminate  is  not  a  mirror  image  of  that 

s 

for  the  [0/90]^  laminate.  Near  the  free  edge,  the  interlaminar 
normal  stress  has  a  maximum  value  at  the  raidplane  of  both  configura¬ 
tions,  i.e.,  at  z/h  =  0,  as  shown  in  Figure  11.  These  results  are 
in  good  agreement  with  those  reported  by  Wang  and  Crossman  [39]. 

In  order  to  permit  comparison  with  other  published  results, 
the  same  problems  were  solved  again  using  the  material  properties 
of  Reference  [24],  which  are  given  below. 


En  =  138  CPa  (20  Msi) 

E22  *  E33  =  14.5  GPa  (2.1  Msi) 

C23  ”  G31  =  G12  =  6  GPa  (0‘85  Msl) 
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Distributions  of  a  and  r  at  z/h  =  1  for  the  [0/90]  and  [90/0] 
z  yz  s  s 

configurations  were  calculated.  For  the  [0/90 ] ^  laminate,  the 

present  results  agreed  well  with  those  reported  by  Renieri  and 

Herakovich  [24].  Both  methods  predicted  that  for  y/b  <  0.5,  the 

classical  laminated  plate  theory  solution  was  approximately 

recovered.  The  maximum  value  of  near  the  free  edge  as  calculated 

with  the  present  three-dimensional  analysis  was  the  same  as  that 

of  Reference  [24].  This  value  was  much  smaller,  however,  than 

that  reported  by  Wang  and  Crossman  [39]  using  the  same  material. 

Basically,  the  results  of  the  present  work  and  those  reported  in 

References  [24,39,88]  were  close  up  to y/b  =  0.5;  but  in  regions  nearer 

the  free  edge,  results  differed.  While  Wang  and  Crossman  [39] 

indicated  a  zero  value  for  a at  the  free  edge,  both  the  present 

method  and  that  of  Renieri  and  Herakovich  [24]  predicted  the  value 

of  a ^  close  to  the  free  edge  to  be  compression.  The  conclusions 

stated  above  for  o  were  found  to  be  true  also  for  the  interlaminar 
z 

shear  stress  distributions.  For  the  [90/0]  configuration,  the 

s 

present  results  again  agreed  well  with  those  of  Reference  [24]  for 
both  components  nf  interlaminar  stress. 


4.3.2  Angle-Ply  Laminates  Under  Mechanical  and  Hygrothermal  Loading 
Angle-ply  laminates  of  the  configuration  [ +0 ] g ,  and  [+0]^, 
where  0  is  the  ply  angle,  were  also  considered  here,  results  being 
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presented  in  Figures  12-14  for  6  =  20°,  30°,  and  70°.  The  following 
observations  apply  to  each  of  these  configurations. 

The  interlaminar  normal  stress  at  midplane  for  the  [+0] 
laminates  increases  from  compression  at  y/b  =  0  to  tension  at 
y/b  =  1,  i.e.,  at  the  free  edge.  For  the  [+0]  laminates, 
decreases  from  tension  at  y/b  =  0  to  compression  at  the  free  edge. 

The  absolute  value  of  at  y/b  =  0  for  the  [+0]  configurations  is 
about  half  that  for  the  [+0]  configurations. 

Results  presented  here  are  in  disagreement  with  those  reported 
by  Pipes  and  Pagano  [26],  who  noted  that  the  presence  of  interlaminar 
stresses  might  be  considered  an  edge  effect  restricted  to  a  small 
region  near  the  free  edge.  In  their  work  of  Reference  [26],  Pipes 
and  Pagano  presented  a  method  for  predicting  the  interlaminar 
stresses  under  uniform  axial  extension.  They  considered  symmetric 
angle-ply  laminates  loaded  by  tractions  applied  only  on  the  ends 
x  =  constant.  A  displacement  field  of  the  form 

u  =  Cx  +  U(y,z) 

v  =  V(y,z)  (56) 

w  =  W(y,z) 

was  then  assumed.  This  type  of  displacement  field,  together  with 
the  symmetry  conditions,  led  to  the  following  boundary  conditions: 


vxz(y,0)  =  0 


(5/) 


Yyz(y.°)  =  0 
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Y  (0 ,  z)  =  0 
yz 

Pipes  and  Pagano  [26]  noted  that  such  an  analysis,  which  considers 
stress  boundary  conditions  only,  cannot  handle  a  bar  with  clamped 
ends  under  extension,  such  as  a  tensile  coupon,  and  that  their 
method  was  simply  an  attempt  to  discern  the  influence  of  a  free  edge 
on  laminate  response.  They  stated  that  the  true  behavior  can  only 
be  studied  by  abandoning  the  assumption  that  all  stress  components 
are  independent  of  x.  Most  of  the  methods  of  analysis  presented  in 
the  literature,  however,  have  taken  the  same  approach  as  Pipes  and 
Pagano  [26]  in  dealing  with  the  problem  of  interlaminar  stresses. 

It  is  believed  that  the  differences  between  the  results  of  Reference 
[26]  and  those  of  the  present  analysis,  which  includes  the  dependence 
on  x,  are  due  to  the  assumptions  used  in  Reference  [26],  as  discussed 
above . 

The  interlaminar  shear  stress  x  ,  presented  in  Figure  13  for 

yz 

different  values  of  0,  is  always  positive  for  the  [+0]  laminates, 

in  agreement  with  results  for  a  [+45 ]  laminate  reported  by  Rybieki 

[34].  The  plots  of  x  shown  in  Figure  13  indicate  a  steady 

v  z 

increase  (or  decrease)  to  a  high  absolute  value  near  the  free  edge 
before  reaching  zero  at  y/b  =  1 ,  as  dictated  by  the  traction-free 
edge  condition.  This  high  value  of  x  near  the  free  edge  was  also 
reported  by  Pipes  and  Pagano  [261.  The  variation  in  the  interlaminar 
shear  stress  x  ,  as  shown  in  Figure  14,  is  less  titan  the  variation 
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Figures  15  through  17  show  the  variation  of  the  interlaminar 
edge  stresses  as  a  function  of  the  lay-up  angle  0  for  a  [+G] 

s 

laminate.  Two  types  of  loading  are  considered,  mechanical  and 
hygrothermal  loadings.  The  present  analysis  can  handle  both 
types  of  loadings,  including  any  combination  of  them.  Hygrothermal 
loading  includes  a  temperature  change,  a  moisture  change,  or  a 
combination  of  the  two.  It  can  be  a  constant  overall  change,  or 
even  some  distribution  of  temperature  and  moisture.  However,  for 
simplicity  and  clarity  in  the  present  discussion,  the  laminate  in 
this  example  will  be  loaded  separately  by  a  mechanical  loading  of 
e"x  =  0.01  percent,  and  a  uniform  thermal  loading  of  AT  -  -50°C,  with 
the  stress  free  temperature  assumed  to  be  room  temperature,  viz, 
21°C.  Apart  from  a  small  variation  due  to  temperature-dependent 
material  properties,  and  a  sign  change,  stresses  produced  by  a 
positive  change  in  temperature  are  the  same  as  those  produced  by  a 
negative  change. 

Several  interesting  results  are  noted  in  comparing  the  free 

edge  effect  of  mechanical  uniaxial  loading  and  thermal  loading,  for 

different  values  of  9.  The  interlaminar  normal  stress  n  ,  shown  in 

z 

Figure  15,  increases  rapidly  as  9  increases  to  about  20°,  and  then 
decreases;  the  decrease  is  much  slower  beyond  0  =  45°.  The  inter¬ 
laminar  shear  stresses  (Figures  16  and  17)  decrease  from  zero  to  a 
maximum  absolute  value  for  values  of  9  in  the  range  of  20°  to  25°. 
Then  they  increase  again  to  zero  at  0  =  90°.  Again,  beyond  9  =  45°, 
the  slope  of  the  curve  decreases  sharply.  The  value  of  0  at  which 
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a)  Uniaxial  Strain  ex  =  0.01%  I 
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a  peak  absolute  value  is  reached  by  the  interlaminar  stresses  varies 
under  the  mechanical  loading  and  thermal  loading.  Under  the  thermal 
loading  condition  of  AT  =  -50°C,  this  peak  value  of  0  is  from  5°  to 
10°  higher  than  for  uniaxial  mechanical  loading. 

Although  moisture  changes  have  not  been  considered  in  this 
example,  the  program  is  capable  of  handling  hygrothermal  loadings 
due  to  moisture  and/or  temperature  changes.  Moisture  changes  are 
handled  by  the  analysis  in  the  same  way  as  temperature  changes. 

4.4  Interlaminar  Stresses  Around  Circular  Holes 

The  presence  of  interlaminar  stresses  near  the  free  edges  of 
laminated  plates  was  demonstrated  in  the  previous  two  examples.  A 
more  complicated  and  interesting  type  of  problem  is  that  of  inter¬ 
laminar  stresses  around  cutouts,  e.g.,  holes,  in  laminated 
composites.  Such  problems,  involving  curved  rather  than  straight 
boundaries,  are  generally  more  difficult  to  analyze.  However, 
since  the  present  analysis  is  a  three-dimensional  finite  element 
approach,  it  is  capable  of  handling  any  type  of  boundary  geometry 
with  similar  ease.  The  following  example  illustrates  this  capability. 

The  laminates  included  in  this  example  will  be  analyzed  first 
under  a  uniaxial  strain  (T  =  0.01  percent)  in  order  to  compare 
results  with  results  using  other  methods.  The  multiaxial  loading 
capability  will  then  be  demonstrated  by  applying  biaxial  loading  to 
the  laminates.  Finally,  inelastic  response  will  be  considered. 

Other  methods  in  the  literature  cannot  handle  cases  other  than 


uniaxial  loading,  since  they  are  based  on  the  uniaxial  displacement 
field  first  introduced  by  Pipes  and  Pagano  [26],  as  discussed 
earlier.  Few  earlier  analyses  considered  nonlinear  material 
behavior  [20,24]; inelastic  laminate  response  was  not  considered  in 
any  of  these  methods. 

Four-ply  laminates  of  the  configurations  [0/90]  and  [90/0]  , 

s  s 

containing  circular  holes,  are  considered.  The  overall  dimensions  of 
the  laminate,  as  given  in  Reference  [43],  are: 

Length,  £  =  203  mm 
Width,  w  =  254  mm 
Ply  thickness,  h  =  7.6  mm 
Hole  radius,  R  =  6.25  mm 

The  three-dimensional  grid  used  is  shown  in  Figure  18;  only  one 
quarter  of  the  upper  two  layers  need  be  considered,  because  of 
symmetry.  In  the  first  part  of  this  example,  for  purposes  of 
comparison,  the  material  assumed  is  the  unidirectional  graphite/epoxy 
composite  used  in  Reference  [43],  the  mechanical  properties  of 
which  are: 
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Figure  18.  Finite  Element  Grid  Used  to  Model  a 

Circular  Hole  in  a  Four-Ply  Laminate. 

The  distribution  of  the  interlaminar  normal  stress  o  is  shown 

z 

in  Figure  19  for  the  two  configurations,  viz,  [0/90J  and  [90/0]  , 

s  s 

under  a  uniaxial  strain  =  0.01  percent,  corresponding  to  an 
average  applied  stress  o^  =  7.1  MPa.  Both  laminate  configurations 
show  a  compressive  value  for  at  0  ■  0°,  and  a  tensile  value  at 
9  ■  90®.  Higher  values  at  9  =  0°  and  90°  are  observed  in  the  [90/0J 

s 
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laminate.  Thus,  the  change  in  stacking  sequence  does  not  produce 
mirror  image  distributions.  Also,  in  contradiction  with  results 
reported  in  Reference  [43],  the  interlaminar  normal  stress  near  the 
free  edge  of  the  hole  does  not  change  sign  as  the  stacking  sequence 
is  changed. 

The  analyses  of  composite  laminates  under  uniaxial  loading  are 

helpful  in  understanding  the  complex  behavior  of  interlaminar  stresses. 

In  actual  service,  however,  laminates  are  often  subjected  to  multi- 

axial  states  of  stress.  The  second  part  of  this  example  illustrates 

the  capability  of  the  present  analysis  to  handle  multiaxial  loading 

situations.  Laminates  with  circular  holes,  previously  analyzed  under 

uniaxial  loading,  are  considered  next  under  varying  biaxial  loadings. 

Table  5  shows  the  interlaminar  normal  stress  at  two  different 

locations,  viz,  0  =  0°  and  0  =  90°,  under  different  biaxial  loading 

conditions.  This  table  has  been  generated  by  varying  the  ratio  of 

T^/T  .  For  the  [90/0]g  laminate,  is  always  positive  at  0  =  90°, 

with  almost  the  same  value.  This  value,  shown  in  column  four  in 

Table  5,  is  reached  rapidly  from  zero  under  F  IT  =0  and  then  attains 

x  y 

a  constant  value  at  about  T  /T  =  ].  Plate  dimensions  and  Poisson's 

x  y 

ratios  appear  to  affect  this  behavior;  further  investigations  3re 

needed  to  fully  understand  this  response.  At  0  =  0°,  o  changes  from 

tension  to  compression  as  T  is  increased  relative  to  c"  .  For  the 

x  y 

[ 0/90 1 ^  configuration,  the  interlaminar  normal  stress  is  always  posi¬ 
tive  at  6  *  0°,  but  the  value  decreases  as  the  ratio  c  /T  increases. 

x  y 

At  0  “  0°,  o  increases  from  compression  to  tension  as  TT  /'c  increases, 
z  x  y 
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In  spite  of  the  fact  that  composite  materials  often  exhibit  large 
amounts  of  inelastic  deformation,  no  method  of  analysis  considered 
the  inelastic  behavior  of  composite  laminates.  To  show  hoe  the  pre¬ 
sent  analysis  can  handle  inelastic  material  behavior,  a  different 
material  system  is  used  in  the  next  part  of  this  example.  The 
material  is  Hercules  AS/3501-6  graphite/epoxy,  the  mechanical 
properties  of  which  were  given  in  Table  3.  The  method  of  analysis 
of  Reference  [43]  did  not  consider  the  inelastic  behavior  and  hence 
no  post  elastic  properties  were  presented. 

The  interlaminar  stresses  for  the  two  cross-ply  configurations 

were  calculated  under  a  biaxial  loading  ratio  e^/ey  =  1.25  (see 

Table  5).  Figure  20  is  a  plot  of  the  interlaminar  stresses  at  the 

midplane,  i.e.,  at  z  =  0,  while  the  interlaminar  stresses  at  the 

interface  between  the  90°  and  0°  plies,  i.e,,  at  z  =  h,  are  shown  in 

Figure  21.  For  both  configurations,  the  interlaminar  normal  stress 

a,  is  dominant  at  0  =  90°,  with  a  higher  tensile  value  in  the  [90/0] _ 

laminate.  The  variation  of  oz  between  0=0°  and  9  =  90°  is  much 

greater  at  the  midplane  than  at  z  *  h.  The  interlaminar  shear  stress 

T„  at  z  *  0  increases  in  absolute  value  to  a  maximum  at  0  *  90°  for 
yz 

both  configurations.  However,  at  z  =  h,  a  peak  absolute  value  is 
attained  at  0  *  45°.  The  interlaminar  shear  stress  tzx  behaves  in  a 
similar  manner  at  z  =  h,  except  for  a  change  in  sign,  and  also  attains 
a  peak  value  at  0  *  45°.  However,  at  z  «  0,  tzx  decreases  in  absolute 
value  towards  a  minimum  at  0  *  90°. 


Figure  20.  Interlaminar  Stresses  Around  a  Circular  Hole  at  the  Midplane 
z  =  0,  in  [0/90 ] s  and  [ 90/0] s  Laminates  Subjected  to  Biaxial 
Loading  ~£  IT  =1.25. 
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Figure  21.  Interlaminar  Stresses  Around  a  Circular  Hole  at  z  =  h,  in 
[0/90]s  and  f 90/0 3 s  Laminates  Subjected  to  Biaxial  Loading 
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To  study  the  inelastic  behavior  of  cross-ply  laminates,  the 

biaxial  loading  was  increased,  keeping  the  ratio  constant.  In  the 

case  of  the  [0/90]  laminate,  yielding  occured  first  at  the  free 
s 

edge  of  the  hole  in  the  inner  90°  ply  at  an  angle  slightly  less  than 
45°,  and  in  the  outer  0°  ply  at  an  angle  slightly  greater  than  45°, 
as  indicated  in  Figure  22.  As  the  load  increased,  the  yield  zone 
moved  across  the  45°  line  in  both  plies,  always  beginning  in  the  inner 
ply,  i.e.,  the  90°  ply.  Yielding  in  the  [90/0]^  laminate  started 
at  the  same  locations  and  at  the  same  load  level.  The  progressive 
growth  of  the  yield  zone  in  each  lamina  for  both  configurations  is 
shown  in  Figure  22.  The  pattern  in  which  the  plastic  zone  propagated 
in  each  ply  did  not  change  as  the  stacking  sequence  changed.  First 
failure  was  found  to  occur  at  the  45°  position,  as  shown  in  Figure 
23.  For  both  configurations,  this  first  failure  was  at  the  free 
edge  of  the  hole,  at  the  miaplane.  The  state  of  stress  around  the 
free  edge  of  cutouts  is  indeed  very  complicated.  A  method  of  analysis 
such  as  the  one  presented  here  is  mandatory  when  designing  laminates 
that  will  not  delaminate  at  free  edges,  whether  at  a  straight  free 
edge  or  around  a  cutout.  The  ability  of  this  method  to  handle 
hygro thermal  loadings  makes  it  of  special  value  in  dealing  with 
polymeric  composites,  which  are  specially  susceptable  to  changes  in 
temperature  and  moisture  levels.  Such  changes  in  environmental 
conditions  often  lead  to  triaxlal  states  of  stress,  which  can  only  be 
handled  by  a  full  three-dimensional  analysis. 
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a)  Propagation  of  yield  zone  in  0°  ply 


b)  Propagation  of  yield  zone  in  90°  ply 
Figure  22.  Propagation  of  Yield  Zone  under  Biaxial  Loading 


Section  5 


SUMMARY  AND  CONCLUSIONS 

The  present  study  has  been  concerned  with  the  three-dimensional 
inelastic  problem  of  generally  orthotropic  laminated  composites, 
including  hygrothermal  effects. 

At  the  present  time  there  are  not  enough  experimental  data  to 
permit  the  construction  of  a  'general*  theory  of  plasticity  for 
composite  materials.  However,  analytical  methods  such  as  the  one 
presented  here  should  continue  to  be  investigated  and  correlated  to 
experimental  data,  in  order  to  direct  future  experimental  programs 
toward  attaining  the  goal  of  a  general  theory. 

The  example  problems  presented  in  the  previous  chapter  demon¬ 
strate  only  a  few  of  the  many  potential  applications  of  the  analysis. 
It  would  be  possible  to  combine  several  of  these  features,  such  as 
multiaxial  loading,  mechanical  and  hygrothermal  loading,  etc.,  to 
model  many  different  situations.  The  purpose  here,  however,  is  to 
give  an  indication  of  the  possible  uses  and  applications  of  the 
program.  Experimental  verification  of  these  capabilities  remains  to 
be  performed. 

The  analysis  is  capable  of  modeling  any  laminated  composite 
subjected  to  triaxial  mechanical  and/or  hygrothermal  loadings.  It 
can  also  be  used  as  a  three-dimensional  micromechanics  analysis, 
since  it  can  model  any  number  of  materials  simultaneously.  Being 
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a  three-dimensional  analysis,  it  is  the  only  way  to  analyze  materials 
such  as  3-D  carbon-carbon  composites,  which  consist  of  yarn  bundles 
oriented  in  three-dimensions.  A  three-dimensional  finite  element 
analysis  is  used  to  model  the  laminate.  The  constituent  material 
properties,  which  can  be  functions  of  temperature  and  moisture 
content,  are  input  as  coefficients  of  a  second  order  polynomial. 

This,  however,  can  be  changed  easily  to  describe  any  other  functional 
relationship  between  the  material  properties  and  temperature  and 
moisture.  The  current  version  of  the  computer  program  uses  a  finite 
element  analysis  based  upon  a  displacement  formulation,  and  linear 
isoparametric  brick  elements  with  three  degrees  of  freedom  per  node. 
Other  types  of  elements  can  be  added  to  the  program  to  model 
geometrically  complex  boundaries,  crack  elements,  or  any  specialized 
elements  such  as  those  of  Reference  [88].  The  modular  form  in  which 
this  program  is  built  makes  it  a  relatively  easy  task  to  do  so. 

Nonlinear  (elastoplast ic)  material  behavior  is  included  by 
means  of  the  tangent  modulus  method.  The  onset  of  plastic  deforma¬ 
tion  is  determined  by  the  use  of  a  yield  surface  which  is  dependent 
on  both  temperature  and  moisture  content  in  the  composite.  Work 
hardening  of  the  material  is  assumed  to  be  isotropic  in  the  present 
version  of  the  program,  for  lack  of  experimental  data  to  describe 
this  phenomenon  otherwise  for  composite  materials.  The  Bauschinger 
effect  is  also  neglected  for  similar  reasons.  Statistical  variations 
in  the  properties  of  typical  composite  materials,  however,  may  well 
prove  to  exceed  the  small  inaccuracies  induced  by  these  assumptions. 
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The  failure  surface  incorporated  in  the  present  analysis  is  similar 
in  shape  to  the  yield  surface  and  is  also  dependent  on  temperature 
and  moisture.  The  failure  of  any  lamina  will  cause  the  current 
computer  program  to  halt  execution,  with  proper  messages.  This  does 
not  always  mean  that  the  laminate  has  failed,  however,  as  discussed 
earlier.  Many  composites  will  continue  to  carry  load  long  after  the 
matrix  has  experienced  many  local  failures  (microcracks) .  For  this 
reason,  crack  propagation  should  be  included  in  future  work.  This 
capability  can  be  incorporated  into  the  present  analysis,  e.g.,  by 
reducing  the  material  stiffness  properties  of  failed  elements. 

Other  methods,  originally  developed  for  unidirectional  composites, 
such  as  those  discussed  in  References  [89-91]  could  also  be  modified 
to  suit  the  laminate  analysis. 

The  addition  of  time-dependent  deformation  into  the  present 
analysis  is  another  feasible  capability  that  could  be  handled. 
Recently  Schaffer  and  Adams  [92]  presented  a  nonlinear  viscoelastic 
t  talysis  for  unidirectional  composites.  In  their  work,  nonlinear 
viscoelastic  constitutive  equations  for  an  isotropic  material  were 
developed  and  incorporated  in  an  elastoplastic  computer  program  [18] 
originally  developed  for  time-independent  analyses.  To  similarly 
add  time-dependent  material  response  to  the  present  work,  visco¬ 
elastic  constitutive  equations  for  a  generally  orthotropic  material 
will  have  to  be  developed. 
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In  summary,  the  present  analysis  is  a  fully  three-dimensional, 
elastoplastic  analysis.  It  can  be  used  to  analyze  any  composite 
laminate  or  structure  consisting  of  any  number  of  generally  ortho¬ 
tropic  materials,  for  which  mechanical  properties  can  vary  with  both 
temperature  and  moisture.  Hygrothermal  and/or  mechanical  loadings 
can  be  handled. 
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APPENDIX  A 


TRANSFORMATION  EQUATIONS  FOR  STRESS 
COMPONENTS  FOR  A  FIBROUS  LAMINA 

Assume  that  the  stress  components  of  a  certain  lamina  are  known 
in  the  reference  axes  (x,y,z),  and  that  the  axes  of  anisotropy  for 
this  lamina  are  (1,2,3)  such  that  the  3-axis  and  z-axis  coincide. 
Then,  if  the  1-axis  is  inclined  at  a  counterclockwise  angle  0  to 
the  x-axis,  the  stress  components  in  the  (1,2,3)  system  can  be  found 
according  to  the  equations 


APPENDIX  B 


FORMULATION  OF  STIFFNESS  MATRIX  FOR  AN 
ISOPARAMETRIC  ELEMENT 

Consider  the  local  coordinates  £,  n,  and  with  the  origin 
located  at  the  center  of  the  element  as  shown  in  Figure  B1 .  These 
coordinates  will  be  so  defined  as  to  give  values  of  +1  and  -1  on  the 
faces  of  the  element.  Thus,  in  this  coordinate  system  the  element 
is  a  cube.  Figure  B2.  In  general,  the  coordinates  of  any  point  are 
defined  by  the  expressions 

8 

x  =  NjXj  +  N2x2  +  ...  +  NgXg  =  l  N^ 

8 

y  =  l  N.y.  (B-l) 

1 

8 

z  -l  tVi 

lii 

where  x^,  y^,  and  are  the  nodal  coordinates  in  (x,y,z)  system  and 
N^  are  the  shape  functions.  Displacements  are  defined  in  the  sane 
way;  hence  the  name  isoparametric. 
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u  =  I  Mu 
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l  NIVi 


v  = 


(B-2) 
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8 

w  -  l  Niwt 
1 

where  u^,  v^,  and  are  nodal  displacements  in  the  x,  y,  and  z 
directions,  respectively,  and  i  ranges  from  1  to  the  total  number  of 
nodes  (8  in  this  example  of  an  eight-node  quadrilateral  element). 

The  shape  functions  are  given  by 

Nt  =  |(1  +  Cti)(l  +  nni)(l  +  CCi)  (B-3) 

where  £^,  n^,  and  ate  the  nodal  values,  i.e.,  +1. 

The  element  matrices  are  formulated  next,  in  terms  of  the 
isoparametric  coordinates.  Relations  between  derivatives  in  the  two 
coordinate  systems  are  established  by  the  chain  rule  of  differentia- 
t  ion , 


(B-4) 


where  commas  denote  partial  differentiation,  and  [J]  is  the  so-called 
Jacobian  matrix.  From  Eq.  (B-l), 
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But  the  strain-displacement  relation  can  be  written  as 


(e  } 


'yz 


V 

l/xy 


u. 


10000000  0: 
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010100000 


V, 


v,y  > 


V, 


W, 


(B-7) 


Substituting  from  Eq.  (B-6)  Into  Eq.  (B-7),  and  utilizing  Eq.  (B-2) , 
the  relation  between  the  strain  and  displacement  is  obtained 


(c)  =  [B] (d } 


(B-8) 


where 


{d )  =/ 


(B-9) 


k  W8  ) 

is  the  displacement  vector. 
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The  element  stiffness  matrix  is  then  formed  as 


[k] 


[BT] [C] [B]d(Vol.) 


(B-10) 


Vol. 


where  [C]  is  the  element  stiffness  matrix  as  given  in  Appendix  C. 


APPENDIX  C 


GENERALLY  ORTHOTROPIC  STIFFNESS  COEFFICIENTS 

The  C  coefficients  in  the  generally  orthotropic  stiffness  matri 


in  Eq.(B- 

10) 

are  listed  below: 
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where 


1/8  (3CU  +  3CU  +  4Cm) 
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APPENDIX  D 


INPUT  INSTRUCTIONS 
FOR  NCLAP 


Card  1  (80A1) 

TITLE  Problem  title 


Card  2  (1115) 

NTP 

NTE 

NPE 

NMAT 

NLI 

NFIX 

NGAUS 

ICHK 

ISTR 


IDIS 


Total  number  of  nodes 
Total  number  of  elements 
Number  of  nodes  per  element 
Number  of  materials 
Number  of  load  increments 
Number  of  constrained  nodes 
Number  of  integration  points 
-  1  Check  input  data 

=  0  No  check 

=  0  Print  element  average  stresses 

=  1  Print  element  stresses  at  every  Gaus  point 

=  2  Full  print  out  of  stresses 

=  0  No  print  out  of  displacements 

=  1  Print  accumulated  displacements 

=  2  Print  incremental  and  accumulated 

displacements 
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THYG 


Cards  3 
NF, 

NOE 
I  MAT 
ANGLE 

IFLOD 

Cards  4 
NP 

COORD 

STEMP 

SMOIS 

Cards  5 

NFIXP 

CODE 

GDIS 


=  1  Print  nodal  temperatures  and  moisture 
distribution 

=0  No  hygrothermal  information  printed 

(1015, F10. 0,15) 

Element  number 

Vector  of  element  connectivity 
Element  material  number 

Element  material  angle,  i.e.,  angle  between 
principal  direction  1  of  the  material  and  the 
x-axis  measured  clockwise 

=  1  Mechanical  load  applied  to  element 
=0  No  mechanical  load  applied 

(15, 5F10.0) 

Node  number 

Nodal  coordinates  (x,y,z) 

Nodal  temperature 
Nodal  moisture  content 


(15, 2 X, 311 , 3F10 . 0) 

Constrained  node  number 
Vector  of  constraining  code 
=  1  Contrained  displacement 

=0  No  constraint 

Prescribed  displacements  in  x,  y,  and  z  directions 


H 


Card  6 
MAT 

Cards  7 
PROPS 

Card  8 
IMLD 

IHCR 


Cards  9 

DTEMP 

DMOIS 

Card  10 
NE 

LNODS 

Cards  11 
LN 


(15) 

Material  number 

(6E10.4) 

Material  properties  (see  Tables  3,4) 

(212) 

Mechanical  loading  flag 

=  1  Increment  contains  applied  mechanical  load 
=0  No  mechanical  load 

=0  No  hygrothermal  load 

-  1  Uniform  change  in  hygrothermal  loading 

=  -1  Hygrothermal  distribution 

(2F10.0) 

Uniform  applied  temperature  change 
Uniform  applied  moisture  change 

(215) 

Mechanically  loaded  element  number 
Number  of  nodes  loaded 

(15, 3F10.0) 

Loaded  node  number 


XL , YL , ZL 


x-load,  y-load,  and  z-load 


Note 
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If  IHGR  on  Card  8  Is  equal  to  -1 .Cards  9  would  then  give 
the  change  of  temperature  and  moisture  for  all  nodal 
points  ending  with  the  last  node.  Cards  9  are  omitted 
totally  if  IHGR  is  zero. 


